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ABSTRACT 


X  , 

The  Doppler  spectrum  of  radar  return  from  the  sea  at  HF  contains  two 
narrow  lines  displaced  upward  and  downward  from  the  carrier  frequency, 
resulting  from  backscatter  off  ocean  surface  waves  moving  toward  and  away 
from  the  radar,  respectively.  These  Bragg  lines  indicate  resonant  back¬ 
scatter  which  occurs  for  the  ocean  waves  of  length  one-half  the  radio 
wavelength.  The  phase  velocity  of  these  gravity  waves  consists  of  two 
components;  the  first  is  determined  by  the  wavelength  and  the  second  by 
the  current  component  in  the  direction  of  wave  propagation  averaged  over 
depth  with  an  exponential  weighting  function  that  has  a  characteristic 
scale  proportional  to  the  wavelength.  The  Doppler  shift  of  the  radar 
carrier  is  determined  by  the  wave  phase  velocity.  Its  variation  with 
carrier  frequency  is  thus  related  to  the  vertical  profile  of  the  current 
component  in  the  radar  direction  by  a  Laplace  transform;  therein  lies 
the  principle  of  radio  measurement  of  ocean  current  and  its  vertical 
shear. 

Radio  backscatter  experiments  to  verify  the  feasibility  of  such 

measurements  were  conducted  at— P^&oadero— oft-the.J2aJ.if or nia  coast,  using 

£. 

the  Stanford^  Radar  operating  at  four  frequencies  covering  the  range  from 
3  to  30  MHz.  The  depth-averaged  radial  current  deduced  from  the  centroid 
of  the  Doppler-shifted  sea  echo  itrt  "the . January  199-8^  experiments  at  6.8, 

13.3,  21.7,  and  29.8  MHz  showed  fluctuations  on  the  order  of  1  cm/sec 

.■>  be  j/ 

superimposed  on  temporal  trends  that  reached  maximum  values  of  x~40  cm/ 
sec.^  Concurrent  in-situ  measurements  were  obtained  by  tracking  spar 
buoys  1,  3,  6,  and  12  m  long,  using  a  microwave  ranging  system.  Both 
radar-  and  buoy-inferred  currents  agree  generally  in  their  variation  with 
time  and  depth;  the  discrepancy  between  them  was  no  larger  than  the  dif¬ 
ference  in  the  currents  measured  by  buoys  deployed  at  different  locations 
within  the  same  range  cell,  typically  from  a  few  to  10  cm/sec.  Inversion 
of  the  Laplace  transform  relating  the  radar-inferred  Doppler  velocity  as 
a  function  of  frequency  to  the  vertical  profile  of  the  current  is  achieved 
numerically  by  seeking  a  solution  that  simultaneously  minimizes  a  mean- 
square  error  and  a  mean-square  deviation  from  an  initial  estimate.  Two 
estimators  of  Bragg-line  width  based  on  the  area  of  the  spectral  curve 
and  the  second  moment  about  its  centroid,  respectively,  are  applied  to 
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The  observed  width  in  units  of  current  is 


found  to  vary  in  inverse  proportion  to  frequency  raised  to  a  power  of  ap¬ 
proximately  0.4.  Its  dependence  on  pulse  width,  range,  and  wind  condi¬ 
tions  is  significantly  weaker. 

— It  is  concluded  that  multifrequency  backscatter  ground-wave  radar  at 
HF  constitutes  a  powerful  technique  for  mapping  current  and  its  vertical 
profile  in  the  top  few  meters  of  the  ocean.  Measurement  resolution  of  a 
few  centimeters  per  second  is  attainable  witti^-a.  coherent  integration  time 
of  ~15  minutes.  Large  areas  of  the  ocean  surface  can  be  scanned  by  range 
gating  and  beam  forming,  with  grid  resolution  on  the  order  of  kilometers. 
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Chapter  I 


INTRODUCTION 


A.  Motivation 

This  research  investigates  the  feasibility  of  remotely  monitoring 
the  ocean  surface  current  and  its  vertical  distribution  using  a  four-fre¬ 
quency  pulsed  radar  operating  in  the  High-Frequency  (HF)  region  of  3  to 
30  MHz.  Both  the  theoretical  and  experimental  aspects  of  this  remote¬ 
sensing  technique  will  be  studied  in  detail. 

The  term  "current"  as  applied  in  this  work  refers  to  the  bulk  motion 
of  the  water  mass  in  the  ocean.  This  motion  is  generated,  in  part,  by 
pressure  and  temperature  differences  and  by  tidal  forces.  It  can  also  t>e 
caused  by  local  wind  blowing  over  the  sea  surface,  with  the  direct  effect 
of  a  wind-drag-generated  shear  current  and  the  indirect  effect  of  a  wave- 
drag-generated  Stokes  current  [Wu,  1975] .  Shear  current  is  expected  to 
follow  a  logarithmic  variation  with  depth,  similar  to  the  logarithmic  wind 
profile  immediately  above  the  sea  surface  [Bye,  1965;  Wu,  1975].  Stokes 
current  is  the  result  of  nonlinearities  in  the  dynamics  of  the  wave,  and 
it  decreases  exponentially  with  depth  [Phillips,  1966]. 

Measurements  of  current  in  the  upper  layers  of  the  ocean  are  impor¬ 
tant  in  the  study  of  momentum  and  mass  transfer  across  the  air/sea  inter¬ 
face,  and  these  processes  significantly  influence  the  global  weather  pat¬ 
tern.  Such  measurements  also  have  direct  practical  applications.  For 
example,  to  minimize  the  adverse  effects  of  oil  spills,  the  direction  and 
magnitude  of  the  drift  of  these  pollutants  must  be  known.  This  drift  is 
partly  the  result  of  the  underlying  ocean  current.  Actually,  any  particle 
that  floats  in  the  ocean  waters  will  be  advected  by  the  bulk  motion  of 
the  water  masses;  some  important  examples  being  fish  larvae  and  food. 

Conventional  methods  of  current  measurement  rely  on  either  moored 
current  meters  or  drifters.  In  the  first  measurement,  the  fluid  velocity 
at  one  fixed  location  is  averaged  over  time  which  results  in  a  Eulerian 
estimate.  The  second  measurement  is  Lagrangian  in  nature  because  each 
drifter  follows  the  mean  drift  of  a  group  of  water  particles.  When  the 
current  field  is  not  uniform  over  the  ocean  surface,  a  large  number  of 
drifters  or  current  meters  must  be  deployed  to  obtain  statistically 


significant  measurements.  The  hostile  sea  environment  renders  such  oper¬ 
ations  extremely  difficult. 

In  the  radio  method,  a  radar  is  employed,  covering  distances  up  to 
40  km.  The  ocean  surface  waves  act  c.s  tracers,  and  their  advection  by 
the  underlying  current  is  detectable  by  the  radar.  The  radar-inferred 
current  is  thus  Lagrangian  in  nature.  Operating  with  narrow  pulses  and 
a  small  beam  width,  the  spatial  resolution  of  a  given  ocean-surface  patch 
illuminated  by  the  radar  is  approximately  10  x  10  km.  By  range  gating, 
portions  of  the  ocean  surface  at  various  radial  distances  from  the  radar 
can  be  examined  simultaneously.  Concurrent  coverage  in  the  azimuthal 
direction  can  be  achieved  by  electronically  steering  the  antenna  beam 
formed  by  an  antenna  array.  After  being  proven  feasible,  this  radio¬ 
measurement  technique  can  be  of  great  utility  to  oceanographers. 


B.  Historical  Background 

Sea  state  here  refers  to  ocean  surface  roughness.  The  ocean  wave 
slope  at  a  given  wavelength  is  commonly  used  as  a  perturbation  parameter 
in  the  analysis  of  both  surface-wave  dynamics  and  the  scattering  of  radio 
waves  from  the  sea  surface. 

Since  Crombie  [1955]  first  identified  the  physical  mechanism  respon¬ 
sible  for  the  first-order  HF  sea  echo,  significant  advances  have  been  made 
in  the  use  of  HF  radar  for  the  remote  sensing  of  sea  state.  A  closed-form 
relation  between  the  radar  cross  section  and  ocean  wave -height  spectrum 
has  been  developed  theoretically  on  the  basis  of  perturbation  analysis 
[see,  for  example,  Barrick,  1970;  Barrick,  1972;  Johnstone,  1975],  and 
the  first-order  results  have  been  experimentally  verified  [Barrick  et  al, 
1974;  Teague  et  al,  1975].  HF  measurements  of  the  wave-height  directional 
spectrum  by  a  synthetic-aperture  technique  [Tyler  et  al,  1974]  and  by  in¬ 
version  of  the  second-order  Doppler  spectrum  [Barrick,  1977a;  Lipa,  1978] 
compared  favorably  with  m-situ  buoy  measurements.  Mapping  of  wind  fields 
in  the  ocean  at  long  ranges  (a  few  thousand  kilometers)  by  means  of  the 
first-order  Doppler  spectrum  has  also  proven  possible  [Long  et  al,  1973; 
Barnum  et  al,  1977]. 

Use  of  the  first-order  sea  echo  to  probe  the  ocean-surface  current 
depends  on  the  fact  that,  in  the  absence  of  current,  the  echo  Doppler 


spectrum  consists  of  narrow  lines  symmetrically  displaced  from  the  radar 
carrier  frequency  by  a  known  amount.  The  observed  small  displacement  of 
the  Bragg  lines  from  the  known  position  is  attributed  to  current  and  is 
referred  to  as  current-induced  Doppler  frequency.  Stewart  and  Joy  [197-?] 
derived  an  integral  relation  between  this  measured  Doppler  frequency  at 
any  radar  frequency  and  the  actual  ocean-current  distribution  with  depth. 
By  simplifying  the  result  to  that  for  a  linear  current  profile,  they  ob¬ 
served  that  the  measured  current-induced  Doppler  frequency  was  equivalent 
to  that  produced  by  uniform  ocean  bulk  motion  at  a  velocity  equal  to  the 
actual  current  at  a  depth  of  approximately  1/24  of  the  radio  wavelength. 
Based  on  this  result,  they  obtained  favorable  comparisons  between  radar- 
inferred  ocean  current  and  in-situ  measurements  by  a  parachute  drogue. 

Barrick  et  al  [1976,1977b]  described  a  current-mapping  HF  radar  sys¬ 
tem  that  measures  the  depth-averaged  ocean  current  at  one  radar  frequency. 
They  extended  a  novel  approach  of  Crombie  [1972]  to  measure  the  direction 
of  signal  arrival  by  comparing  the  phases  of  the  signal  on  adjacent  re¬ 
ceiving  antennas.  This  method  eliminates  the  need  for  a  large  antenna 
array  to  obtain  spatial  resolution  in  the  azimuthal  direction.  They  have 
been  able  to  produce  two-dimensional  maps  of  radar-inferred  ocean  current 
over  a  portion  of  the  Gulf  Scream  [Barrick  et  al,  1977b]  and  at  the  lower 
Cook  Inlet  in  Alaska  [Barrick,  1978a] . 

C.  Contributions 

The  major  contributions  of  this  investigation  are  as  follows: 


•  critical  examination  of  the  perturbation  approach  used  by 
Stewart  and  Joy  [1974]  to  derive  the  ocean-wave  dispersion 
relationship  in  the  presence  of  a  current  that  varies  with 
depth;  their  first-order  result  is  extended  here  to  second 
order  for  an  exponential  drift  profile 

•  derivation  of  the  first-order  Bragg-scatter  theory  for  a 
pulsed  monostatic  radar,  based  on  simple  wave-propagation 
concepts,  to  demonstrate  the  current-measuring  capability 
of  backscatter  radar  and  to  investigate  the  limitations  of 
the  radar  technique  in  a  straightforward  manner 

•  analysis  of  data  collected  at  Pescadero  on  the  California 
coast  in  January  1978,  using  a  four-frequency  HF  backscatter 


3 


radar  to  produce  ocean-current  estimates  with  a  consistency 
of  a  few  centimeters  per  second,  corresponding  to  an  accu¬ 
racy  of  typically  a  few  percent 


•  deduction  of  the  vertical  profile  of  ocean  current  from  the 
radar  data  at  four  frequencies  and  comparison  to  the  profile 
estimated  from  in-situ  drift  measurements  of  spar  buoys  of 
four  different  lengths,  assuming  that  the  profile  is  loga¬ 
rithmic  with  depth;  similar  features  in  the  deduced  profiles 
have  been  observed 

•  development  of  a  stabilized  numerical  inversion  algorithm 
to  obtain  the  vertical  profile  of  ocean  current  from  the 
Doppler  velocity  measurements  at  four  frequencies  without 
requiring  the  profile  to  be  logarithmic  in  depth 

•  study  of  the  behavior  of  the  width  of  the  Bragg  lines  based 
on  over  100  hours  of  data  collected  at  Pescadero,  using  a 
four-frequency  HF  backscatter  radar,  spanning  a  period  of 
four  years 


D.  Organization 


Chapter  II  investigates  the  hydrodynamics  of  waves  on  an  ocean  cur¬ 
rent  that  varies  with  depth.  The  ocean-wave  dispersion  relation  is  de¬ 
rived,  based  on  basic  principles  and  the  perturbation  approach  of  Stewart 
and  Joy  [1974] .  Their  perturbation  result  is  extended  to  second  order 
for  an  exponential  current  profile. 

The  concept  of  ocean-wave  spectra  is  discussed  in  Chapter  III.  Vari¬ 
ous  definitions  of  the  directional  wave  spectrum  are  examined  and  compared 

In  Chapter  IV,  the  first-order  Bragg  scatter  theory  is  derived,  based 
on  simple  wave -propagation  concepts .  The  result  is  then  used  to  consider 
the  effects  of  an  ocean  current  that  exhibits  horizontal  and  vertical  va¬ 
riations. 

Chapter  V  describes  the  Stanford  four-frequency  radar  system  and  the 
signal-processing  procedure  involved  in  producing  Doppler  spectra  from  the 
radar  data  collected.  A  filtering  scheme  is  developed  to  isolate  the 
first-order  Bragg  signal  from  the  second-order  effects  and  system  noise. 
Included  is  a  discussion  of  possible  system-introduced  artifacts  that  are 
shown  to  have  negligible  effects  on  the  data.  In-situ  current  measurements 
obtained  by  drifters  are  also  described. 
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The  centroid  estimator  of  the  current-induced  Doppler  velocity  is 
analyzed  in  Chapter  VI.  It  is  demonstrated  that,  under  appropriate  con¬ 
ditions,  the  estimator  produces  an  unbiased  and  consistent  measure  of  the 
ocean-current  component  in  the  direction  of  the  radar  averaged  over  depth 
and  over  the  radar  range  cell.  The  temporal  characteristics  of  the  aver¬ 
aged  current  thus  estimated  will  be  examined  from  four  days  of  data  col¬ 
lected  at  Pescadero  in  January  1978. 

The  problem  of  recovering  the  actual  current  distribution  with  depth 
from  the  radar -inferred  averaged  current  at  the  four  radar  frequencies  is 
studied  in  Chapter  VII.  A  logarithmic  current  profile  is  assumed  to  fa¬ 
cilitate  comparisons  between  che  radar-inferred  and  buoy-inferred  profiles. 
The  instability  problem  in  the  numerical  inversion  of  the  integral  equa¬ 
tion  needed  to  reconstruct  the  current  profile  from  the  radar-inferred 
depth-aver  2d  current  is  then  examined;  an  algorithm  constructed  to  re¬ 
solve  this  problem  is  computerized,  and  its  performance  is  evaluated. 

Chapter  VIII  presents  the  estimated  bandwidths  of  the  Bragg  lines 
for  the  entire  set  of  radar  data  obtained  from  May  1975  through  January 
1978.  The  behavior  of  the  observed  width  as  a  function  of  radar  frequency 
and  pulsewidth,  range  from  radar,  and  wind  speed  and  direction  is  studied 
to  achieve  an  understanding  of  the  underlying  causes  of  the  finite  Bragg 
width. 

Chapter  IX  is  a  summary  of  the  results  of  this  work,  and  areas  for 
future  research  are  recommended. 


Chapter  II 


PHASE  VELOCITY  OF  FIRST-ORDER  OCEAN  GRAVITY  WAVES  IN 
THE  PRESENCE  OF  OCEAN  CURRENTS  WITH  VERTICAL  SHEAR 

As  will  be  demonstrated  in  Chapter  IV,  HF  radars  can  detect  ocean 
surface  waves  having  wavelengths  from  a  few  meters  to  a  few  tens  of  me¬ 
ters.  The  characteristics  of  these  waves  are  not  significantly  controlled 
by  capillary  forces;  their  natural  modes  are  determined  only  by  gravity 
and,  for  this  reason,  they  are  called  "gravity  waves." 

This  investigation  is  concerned  primarily  with  gravity  waves  whose 
average  height  h  (the  height  from  crest  to  trough)  is  small  compared 
to  their  length  L  (crest-to-crest  distance) ,  especially  those  with 
longer  wavelengths  because  waves  over  a  few  meters  high  are  not  frequently 
encountered.  Ocean  waves  also  break  if  h/L  is  larger  than  approximately 
1/7  [Kinsman,  1965].  The  discussion  in  this  chapter,  therefore,  is  lim¬ 
ited  to  first-order  waves,  and  all  terms  higher  than  first  order  in  h/L 
are  neglected  in  the  mathematical  derivations  to  follow.  Consequently, 
the  resulting  equations  are  linear  in  the  wave-perturbed  quantities  so 
that  the  individual  plane  waves  (called  long-crested  waves)  can  be  con¬ 
sidered.  Only  the  effects  of  perturbations  that  have  a  single  horizontal 
wave -number  vector  k  and  one  radian  frequency  u)  need  be  investigated. 

This  analysis  is  limited  to  deep-water  waves.  Because  the  ocean 
depth  is  assumed  to  be  much  larger  than  the  wavelength  L  of  interest, 
the  ocean  bottom  is  considered  to  be  infinitely  far  away  from  the  air/sea 
interface. 

The  terms  "fluid  element"  and  "water  particle"  are  used  interchange¬ 
ably.  A  fluid  element  is  a  volume  element  small  enough  to  be  considered 
a  particle  and  is  large  enough  to  contain  sufficient  molecules  so  that 
the  concepts  of  local  density  and  local  pressure  are  applicable. 

Wind  is  one  of  the  dominant  forces  that  generates  ocean  waves.  The 
emphasis  of  this  study,  however,  is  not  focused  on  the  wave-generation 
mechanisms  but  on  the  natural  modes  of  gravity  waves.  On  the  other  hand, 
understanding  of  the  wind  generation  of  waves  depends  on  the  measurements 
of  relevant  ocean  parameters,  among  which  is  the  wind-produced  current 
shear  in  the  upper  layers  of  the  ocean.  A  mean  horizontal  current  U(z) 
is  assumed  to  exist  in  the  surface  layers — produced  by  wind  drag  or 
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otherwise.  Here,  £  is  the  positive  vertical  axis,  where  z  =  0  rep¬ 
resents  the  mean  sea  surface  (see  Fig.  2.1).  It  will  be  demonstrated  in 
Section  A  that  only  the  component  of  current  in  the  direction  of  wave 
propagation  will  influence  the  wave  phase  velocity. 


HORIZONTAL 
CURRENT  U 


Fig.  2.1.  GEOMETRY  OF  WAVE/CURRENT 
INTERACTIONS. 


In  the  following  calculations,  U  is  assumed  to  be  constant  in  time 
t  and  in  the  horizontal  coordinates  (x,y) .  The  results  obtained  under 
this  assumption  are  valid  for  a  wave  period  short  compared  to  the  actual 
time  scale  of  the  ocean-current  variation  and  for  a  wavelength  small  com¬ 
pared  to  the  scale  length  of  the  ocean-current  variation. 

Viscous  effects  are  neglected;  they  are  important  only  when  the  Rey¬ 
nolds  number  R  is  comparable  to  or  smaller  than  50  (R  =  UL/V,  where 
U  is  a  representative  fluid  velocity,  L  is  a  representative  spatial 
scale,  and  V  is  the  kinematic  viscosity) .  For  water,  R  is  very  large 
compared  to  unity  in  most  flow  systems  of  importance  [Batchelor,  1970] . 

Based  on  the  above  assumptions,  the  expression  for  the  ocean-wave 
phase  velocity  c  is  derived  in  the  following  section  in  terms  of  the 
wave  number  k  and  ocean  current  U.  Beginning  with  the  equation  of 
motion,  the  inviscid  Orr-Sommerfeld  equation  for  the  vertical  component 
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of  the  wave-induced  fluid  velocity  will  be  derived.  Boundary  conditions 
at  the  free  surface  will  be  linearized,  and  then  a  solution  is  obtained 
for  the  case  where  the  vertical  profile  of  the  ocean  current  is  linear. 
The  problem  of  relating  the  wave  phase  velocity  to  a  general  current  pro 
file  is  solved  by  a  perturbation  approach  wherein  the  current  profile  is 
assumed  to  have  a  small  vertical  shear. 


A.  Derivation  of  the  Orr-Sommerfeld  Equation 

The  restoring  forces  on  gravity-wave  systems  are  pressure  gradients 
and  gravitational  pull.  The  equation  of  motion  for  each  water  particle 
is 


d 

—  v 
dt  -T 


VPm 

—  “  V(gz) 


(2.1) 


where 

vT  =  fluid  velocity  in  an  inertial  frame 
=  scalar  pressure  in  the  water 
g  =  gravitational  acceleration 
p  =  water  density 
V  =  vector  gradient 

z  =  vertical  Cartesian  coordinate  of  the  water  particle 

Note  that  these  restoring  forces  are  all  irrotational .  The  equation  of 
motion  yields  the  following  form  when  the  curl  of  both  sides  is  taken: 


Vx 


0 


(2.2) 


Because  v describes  the  spatial  velocity  distribution  in  the  water, 
acceleration  of  the  particle  is 


d 

—  v 
dt  -T 


9t  -T  +  <?T 
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-T 


(2.3) 
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The  second  term  results  from  acceleration  caused  by  the  migration  of 
water  particles  to  a  region  cf  higher  fluid  velocity. 

The  velocity  field  can  be  separated  into  two  additive  parts — a  con¬ 
tribution  U  that  is  time  independent  and  a  contribution  v  that  is 
not;  therefore, 

vm  =  U  +  v 
-T  —  — 

where  U  is  the  time-averaged  fluid  velocity  to  be  called  ocean  current 
or  mean  water  drift.  This  discussion  is  also  limited  to  ocean  currents 
of  such  large  horizontal  scales  that  their  local  variations  in  the  x,y 
directions  are  negligible  and,  as  a  result, 

U(x,y,z)  =  U(z) 


The  remaining  contribution  v  is  the  wave-induced  perturbation  in 
the  fluid  velocity  decomposed  into  plane-wave  components. 


v  =  v  (z)  exp[i(k  •  r  -  cot)] 


(2.4) 


where 


k  =  horizontal  wave-number  vector  of  a  specific  plane  wave 
<o  =  radian  frequency  of  the  plane  wave 

r_  =  horizontal  vector  xx  +  yy,  where  x,y  are  the  unit  vectors  in 
the  x,y  directions,  respectively 

i  = 

Here,  the  averaged  wave  height  is  small  compared  to  wavelength  so  that 
the  terms  in  Eq.  (2.2)  of  a  higher  order  than  v  can  be  neglected.  Be¬ 
cause  this  assumption  linearizes  Eqs.  (2.2)  and  (2.3),  the  plane-wave 
components  can  be  considered  individually. 

In  this  section,  the  x  coordinate  axis  is  chosen  to  lie  parallel 
to  the  direction  of  the  ocean  current  and,  assuming  purely  horizontal 
current , 


—  v  =  -i  (co  -  k  •  U)  V  +  xw  -r—  U 
dt  -T  —  —  —  dz 
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(2.5) 


II 

«£,v 
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where  u,v,w  are  the  three  Cartesian  components  of  v.  After  some  alge¬ 
braic  manipulation,  the  linearized  version  of  Eq.  (2.2)  is 


(0)  -  k  •  U)  Vxv  =  ( —xv  +  yu)  (k_  •  U ' )  -  iy(wU"  +  w'U')  -  z  (k  •  y)  wU' 

(2.6) 

where  the  primes  indicate  differentiation  with  respect  to  the  vertical 
coordinate  z.  This  vector  equation  contains  three  scalar  components, 
and  only  two  are  independent.  For  convenience,  the  x  and  y  compo¬ 
nents  are  chosen, 


(to  -  k  U)  fik  w  -  v’)  +  k  U'v  =  0  (2.7a) 

x  y  x 

(to  -  k  U)  (u'  -  ik  w)  -  k  U'u  +  i(wU"  +  w’U’)  =  0  (2.7b) 


wheie  k^,k^  ara  the  components  of  k  in  the  x,y  directions,  respec¬ 
tively. 

The  continuity  condition  is  that  the  rate  of  change  of  density  within 
an  infinitesimal  volume  must  ke  balanced  by  the  net  incoming  mass  flow, 
and  this  condition  is  expressed  as 


9p 

3t 


-V 


(PYt> 


which  can  be  rewritten  as 


dp 

dt 


+ 


V 


-T 


U 


where  dp/dt  is  the  rate  of  change  of  density  of  a  fluid  element  moving 
with  the  water  particles.  Because  the  fluid  is  essentially  incompressible 
under  motions  induced  by  gravity  waves  of  decameter  wavelength,  this  rate 
of  change  is  zero  which  results  in  the  incompressibility  condition. 


V  •  vm  =  V  •  v  =  0 

-T  — 

The  second  part  of  this  equation  is  a  result  of  the  assumption  that  the 
current  U  is  horizontal  and  varies  with  depth  only.  For  each  plane-wave 
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component  with  wave  number  k  =  (k^,k^),  the  above  incompressibility 
condition  can  be  expressed  as 


-ik  u  =  ik  v  +  w’ 
x  y 


Note  that  this  is  the  continuity  equation  for  incompressible  fluid,  used 
here  to  eliminate  U  and  U’  from  Eqs.  (2.7).  After  straightforward 
manipulations , 


k  U'v 

ik  w  -  v'  +  —  =  0 

y  A 


ik  (k  U'v)  k  U"w 

ik  v'  +  w"  -  k  w  -  — ~ —  =  0 
y  x  A  A 


(2.9a) 


(2.9b) 


where  A  =  (0  -  k  U . 

x 

Equation  (2.9a)  is  now  multiplied  by  ik^  and  the  result  is  added 
to  Eq.  (2.9b).  The  resulting  second-order  differential  equation  for  the 
vertical  component  of  the  wave-induced  velocity  perturbation  is 


/  k  u"  A 

w"  +  (  -  ■  ■  -  k'J  »  =  0 

\“ '  V  / 


Because  U  and  the  x  axis  are  parallel. 


k  U  =  k  •  U 
x  —  — 


arid  Eq.  (2.10)  can  be  rewritten  as 


(2.10) 


w"  +  - 


(k  •  U)  " 

(U)  -  k  •  U) 


-  k  I  w  =  0 


(2.11) 


which  is  a  form  of  the  Orr-Sommerfeld  equation  for  inviscid  fluid.  Be¬ 
cause  the  scalar  product  is  invariant  under  coordinate  rotation,  it  can 
be  concluded  that,  for  an  arbitrary  ocean-current  direction  in  the  hori¬ 
zontal  plane,  this  equation  is  still  valid.  It  should  be  noted  that  only 


the  component  of  the  mean  horizontal  drift  in  the  direction  of  wave 
propagation  affects  the  coefficients  of  Eq.  (2.11),  which  implies  that 
the  ocean-wave  dispersion  relation  w=w(k)  is  independent  of  currents 
orthogonal  to  the  propagation  direction.  It  is  thus  possible  to  align 
the  wave-number  vector  k  with  one  of  the  horizontal  axes  and  reduce  the 
wave-propagation  problem  to  two  dimensions;  such  a  reduction  is  attrib¬ 
uted  to  Squires  [1933] . 


B.  Boundary  Conditions  at  the  Free  Surface 

The  dynamic  boundary  at  the  ocean  surface  can  be  determined  by  exam¬ 
ining  the  vertical  component  of  Eq.  (2.1).  Based  on  the  assumption  of 
purely  horizontal  mean  drift  with  only  a  vertical  gradient,  this  compo¬ 
nent  simplifies  to 

1  9Pt 

-i  (to  -  k  •  U)  w  = - r - g  (2.12) 

—  —  p  dz  * 


for  each  plane  wave. 

The  pressure  PT  can  be  separated  into  two  components, 


P 


T 


po  +  pi 


where  PQ  is  the  hydrostatic  pressure  in  the  absence  of  waves,  and  P^ 
is  the  dynamic  pressure  induced  by  surface  waves.  The  zeroth-order  com¬ 
ponent  of  Eq.  (2.12),  therefore,  is 


PA(z)  =  P  -  Pgz 

U  ci 


where  z  is  measured  from  the  mean  sea  surface  (positive  upward) ,  and 
P^  is  the  atmospheric  pressure  considered  here  as  an  ambient  condition 
because  no  attempt  is  made  to  analyze  the  generation  of  waves  by  wind 
turbulence . 

The  first-order  component  can  be  obtained  by  multiplying  each  side 
of  Eq.  (2.12)  by  f|  which  is  the  surface  height  above  the  mean  level 
z  =  0  (see  Fig.  2.2)  and  dropping  terms  nonlinear  in  the  wave -perturbed 
quantities.  The  result  is 
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=  Pgn 


z=0 


(2.13) 


which  is  the  dynamic  condition  at  the  free  surface  z  =  n#  where  r)  is 
assumed  to  be  small  enough  for  the  vertical  pressure  gradient  to  be  ap¬ 
proximated  linearly  as 
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Pressure  is  continuous  at  the  interface  z  =  r\.  Because  no  wind/sea  in¬ 
teraction  is  assumed,  there  is  no  first-order  pressure  at  the  surface. 
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Fig.  2.2.  SEA  SURFACE  IN  THE  PLANE  OF 
WAVE  PROPAGATION. 


The  second  condition  at  the  air/sea  interface  is  the  restriction  of 
an  unbroken  surface.  If  z  describes  the  vertical  coordinate  of  a  water 
particle,  this  condition  can  be  stated  [see  Lamb,  1932]  as 


(z  -  n) 

J  z=n 


=  o 


Now,  dz/dt  is  the  vertical  particle  velocity  w.  This  condition  indi¬ 
cates  that  the  vertical  velocity  of  any  particle  at  the  surface  equals 
the  vertical  velocity  of  the  surface  seen  by  the  same  particle, 

w|z=0  =  -i[co  -  k  *  0(0)]  n  (2.14) 
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This  is  called  the  kinematic  boundary  condition  wherein  the  terms  nonlin¬ 
ear  in  the  wave-perturbed  quantities  have  been  dropped.  Combining  the 
dynamic  and  kinematic  conditions  to  eliminate  the  surface-height  variable 
n  results  in 

-i[w  -  k  •  U(0>]  P  =  pgw|z=0 

z=0 

The  x  component  in  Eq.  (2.1)  is  now  used  to  eliminate  the  dynamic  pres¬ 
sure  P^.  After  some  straightforward  algebra, 

-i [to  -  k  ♦  U(0)12  u  |  +  w|  U'  (0)[u>-k  •  0(0)]  =  gk  w|  (2.15) 

For  simplicity,  the  x  axis  is  now  chosen  to  lie  parallel  to  the 
wave-number  vector  k.  The  component  of  ocean  current  in  the  x  direc¬ 
tion  is  denoted  by  U  . 

x 

The  continuity  equation  for  incompressible  fluid  [Eq.  (2.8)]  requires 

that 

iku  =  -w' 

which,  when  substituted  into  (2.15) ,  obtains  the  linearized  boundary  con¬ 
dition  at  the  free  surface, 

|w-kU  (0)1  2  w'  I  +  kU'  (0)  L- kU  (0)  w|  =  gk2w|  A  (2.16) 

L  x  J  'z=0  x  L  x  J  'z=0  J  'z=0 

C.  Wave  Phase  Velocity  Determined  from  Solutions  of  the  Orr-Sommerfeld 
Equation 

The  phase  velocity  c  of  the  ocean  surface  wave  with  wave  number  k 
and  radian  frequency  0)  is 


0) 


As  in  Section  B,  k  is  in  the  positive  x  direction  and  is  the  x 

component  of  the  total  current  vector.  The  Orr-Sommerfeld  equation  [Eq. 
(2.11)]  for  the  vertical  component  w  of  the  wave-induced  particle  mo¬ 
tion  can  now  be  written  as 


w" +  \o~nr  ~ k)  w  -  0 


(2.17) 


where  the  double  primes  indicate  differentiation  with  respect  to  the  dis¬ 
tance  z  above  the  mean  sea  level.  The  linearized  boundary  condition  as 
defined  in  Eq.  (2.16)  is  now 


(c  -  U  J  w'  +  U'  (c  -  U  )  w  =  gw 
xO  0  xO  xO  0  0 


(2.18) 


where  0  denotes  evaluation  at  the  mean  sea  level  z  =  0. 

Equations  (2.17)  and  (2.18)  are  solved  below  for  the  following  cases. 

It  should  be  noted  that,  in  cases  1  and  2,  U"  -  0  so  that  the  solution 

x 

to  the  Orr-Sommerfeld  equation  [f2.17)]  is 


w  =  a  exp(kz) 


(2.19) 


because  w  =  0  at  the  rigid  ocean  bottom  z  =  -°°. 


Case  1:  A  Constant  Current  Profile  U  (z)  =  U„. 

In  this  case,  an  expression  is  obtained  for  wave  phase  velocity  c^ 
in  the  absence  of  current,  and  the  intuitively  obvious  Doppler  shift  in 
the  phase  velocity  in  the  presence  of  a  uniform  current  is  introduced. 

Eased  on  the  solution  in  Eq.  (2.19),  the  linearized  boundary  condi¬ 


tion  becomes 


(c  -  V2  =  i 


(2.20) 


If  there  is  no  mean  drift,  the  gravity  wave  with  wave  number  k  has  a 
phase  velocity  of 


c  =  ±c. 


=  *(k) 


(2.21a) 


where  +  and  -  indicate  waves  propagating  in  the  positive  and  negative  x 
directions,  respectively. 


In  the  presence  of  a  uniform  horizontal  drift  U^,  the  wave  radian 
frequency  is  Doppler  shifted  by  an  amount  equal  to  kU^ — a  result  obtain¬ 
able  from  (2.20), 


c  =  ±c  +  U 
0  0 


(2.21b) 


When  waves  move  with  the  ocean  current,  the  sum  of  c„  and  U  is  ob- 

0  0 

served  and,  when  they  move  opposite  to  the  ocean  current,  the  difference 

of  c„  and  U_  is  observed. 

0  0 


Case  2:  A  Linear  Current  Profile  U  (z)  =  tL  +  (IL/&)  z. 

- - - — ■ — - x  -  0  -  1  - 

This  is  the  only  profile  that  yields  an  exact  expression  for  the  wave 

velocity  c. 

Again,  Eq.  (2.19)  is  the  general  solution  to  the  Orr-Sommerfeld  equa¬ 
tion.  The  boundary  condition  [(2.18)]  reduces  to  the  quadratic  equation. 


<c  -  V2  +  s  <°  -  V  -  lco>2  ■ 0 

whose  solution  is 


c 


(2.22) 


where  s  =  [l/(2k£)]  •  u  /c  =  [oy(2kc  )). 

If  the  term  quadratic  in  U^/cQ  neglected,  the  result  is  identi¬ 
cal  to  that  obtained  by  Stewart  and  Joy  [1974]  who  observed  that  the  dif¬ 
ference  between  c  and  the  still-water  phase  velocity  cQ  is  the  value 
of  the  current  at  depth  l/2k  below  the  mean  sea  level.  If  the  quadratic 
term  is  dropped,  the  magnitude  of  this  difference  is  the  same  for  waves 
propagating  upstream  and  downstream  relative  to  the  current.  Equation 
(2.22)  indicates  that,  for  waves  propagating  downstream,  this  difference 
is 


6 


+ 


1/2 


-  c 


0 


z=-l/(2k) 


+  co(1  +  s2) 


and,  for  waves  propagating  upstream,  it  is 


6  =  U 


-1/ (2k) 


’of1  *  •’) 


(-Cq) 


Waves  propagating  downstream,  therefore,  are  Doppler  shifted  more  than 
those  propagating  upstream  by  an  amount  A  (rad/sec) , 


A  =  k(<S+  -  6_)  =--  2kc 


o  [f1  +  s2) 


(2.23) 


Because  the  quadratic  term  in  Eq.  (2.22)  always  adds  to  c^,  the  waves 
are  speeded  up  by  that  term. 


Case  3:  General  Current  Profile  in  the  Weak- Interaction  Regime. 


In  this  case,  the  maximum  magnitudes  of  U  /c  ,  U'/kc  and  U"/ 

2  x  0  x  0  x 

k  cQ  are  all  on  the  order  of  £  where  0  <  e  «  1.  Following  Stewart 

and  Joy  [1974] ,  the  perturbation  approach  is  adopted.  Their  first-order 

result  wili  be  confirmed  by  the  following  independent  derivations,  and  a 

second-order  result  will  be  obtained  for  an  exponential  profile. 

Equation  (2.17)  will  now  be  solved  by  perturbation  expansion.  Only 

cases  are  considered  where  the  maximum  magnitudes  of  U  ,  U'/k,  and  U"/ 
2  xxx 

k  are  small  compared  to  the  wave  phase  velocity  cQ  in  the  absence  of 

current.  The  perturbation  expansions  for  the  phase  velocity  c  and  for 

the  vertical  component  w  of  the  wave-induced  particle  velocity  are 


(0)  (1)  (2) 
G~C  -f  q  +  c  '  * 


(0)  (1)  (2) 
w  =  w  +  w  +  w  + 


and  the  perturbation  parameter  is 

e  =  max  |r|  0  <  e  <  1 

where  r  =  U  /c_. 


The  n  derivative  of  is  assumed  to  have  a  magnitude  of  order 

kn£o0.  Writing  Pi  =  c(l)/c0,  with  pQ  =  1,  Eqs.  (2.17)  and  (2.18) 
can  be  expressed  as 

w"  -  k2w  =  -r" jl  -  (Pl  -  r)  +  ...j  w  (2.24a) 

with  the  condition  that 


r  -i2 

1 

[1  +  (Px  -  r)  +  p2  +  ...  w*  +  r' 

1  +  (P1  -  r)J 

w  =  kw 


(2.24b) 


at  z  =  0. 

To  first  order  in  £,  Eqs.  (2.24)  simplify  to 


(1)"  ,2  (1) 

w  -  k  w  =  -r"a  exp(kz) 


with  the  condition  that 


r  'D  *  di 

-|WV  -  kw  '  +  (r*  -  2kr)  a 


2ka 


(2.25a) 


(2.25b) 


at  z  -  0.  The  second- order  terms  can  be  similarly  extracted: 


(2)"  2  (2) 
w  -  k  w 


=  r"  -  r)  a  exp(kz)  -  w(1*J 
with  the  condition  at  z  =  0  that 


.it 


(2)'-kw(2) 


(2.26a) 


+  2  (p^  r)  w^  +  r'w^  ^  +  ka(p^-r)2  +  ar'(p^-r)J 


2ka 


(2.26b) 


First-Order  Solutions 


The  general  solution  to  Eq.  (2.25a)  is  the  sum  of  the  complementary 
solution  d  exp(kz)  and  the  particular  solution  w^(z',  where  d  is 


{mw*  W  4,»  «■•  • 


an  arbitrary  constant.  The  solution  proportional  to  exp(-kz)  is  re¬ 
jected  because  it  diverges  at  z  =  -00.  The  particular  solution  can  be 
written  as 


w  =  -b(z)  a  exp(-kz) 
P 


(2.27) 


This  is  substituted  into  Eq.  (2.25a)  to  arrive  at  the  following  differ¬ 
ential  equation  for  b(z): 


whose  solution 


b"  -  2kb'  -  r"  exp(2kz)  =  0 


b'  =  r'  exp(2kz) 


can  be  obtained  by  inspection. 

The  first-order  solution  to  the  Orr-Sommerfeld  equation  is 


(1)  fZ  Ux(Z) 

w  =  d  exp(kz)  -  a  exp(-kz)  I  — -  exp(2kz)  dz 

•'-CO  C0 


The  first  term  will  not  affect  the  first-order  phase  velocity  c.  Direct 
substitution  of  the  above  expression  for  w^  into  the  boundary  condi¬ 
tion  in  Eq.  (2.25a)  yields 


(1)  „ 
c  =  U 


/% 

•'—00 


' (z)  exp(2kz)  dz 
x 


Alternately,  after  integration  by  parts, 


c(1)  -  2k 


r° 

J  u> 

-'—00 


(z)  exp(2kz)  dz 


(2.28) 


which  is  the  same  result  obtained  by  Stewart  and  Joy  [1974]  except  that 
their  final  result  has  misplaced  parentheses  and  their  U  is  the  norma¬ 
lized  current  r. 

It  can  be  concluded  that  the  current-induced  Doppler  velocity  is  the 
ocean-current  component  averaged  over  depth  with  exp(2kz)  as  the 


weighting  function.  The  integral  in  Eq.  (2.28)  can  also  be  interpreted 
as  a  Laplace  transform  evaluated  on  the  real  axis.  These  observations 
will  be  used  to  extract  ocean-current  and  current-shear  information  from 
HF  backscatter  data. 


•  « 


Second-Order  Solutions 


Because  the  general  solution  to  Eq.  (2.25b)  is  complicated  and  dif¬ 
ficult,  only  a  special  case  is  presented  to  illustrate  the  method  and 
provide  an  interesting  result. 

A  current  profile 


U^(z)  =  Uq  exp(mz) 


-1 


decaying  with  depth  and  with  a  scale  length  of  m  is  assumed.  After 
much  algebra. 


2kU  c  „ 

0.02  k 
c  =  c„  +  — — ; +  —  s 


0  2k  +  m  2  m  +  k 


(2.29) 


where  s  =  (U0/cQ)  *m/(2k  +  m),  which  is  valid  for  waves  propagating  in 


the  direction  of  (downstream) .  For  waves  propagating  upstream,  cQ 


is  replaced  by  -cQ,  and  the  phase  velocity  of  these  waves  is 


2k00  °0  2  k 

2  —  ~C  T  - -  ■  "  ’  —  - CJ  - 

'upstream  0  2k  +  m  2  m+k 


To  verify  this  result,  consider  the  case  where  the  magnitude  of  m 
is  much  smaller  than  k;  in  other  words,  the  current  profile  varies 
slowly  (on  the  order  of  one  wavelength)  on  the  spatial  scale.  This  pro¬ 
file  and  the  phase  velocity  can  then  be  approximated  by  a  Taylor  series 
expansion,  thus  yielding 


c  =  U 

z=- 


z=-l/(2k) 


±c0(l+is2) 


(2.30) 


where  s  =  U'  _y(2kc„). 

z=u  0 
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For  a  linear  current  profile,  Eq.  (2.22)  is  an  exact  expression  for 

phase  velocity  c.  If  the  square-root  term  is  expanded,  the  result  will 

2 

agree  with  Eq.  (2.30)  to  order  s  . 


D.  Summary 


Still-water  phase  velocity  of  ocean  surface  waves  of  decameter  wave¬ 
length  has  been  shown  to  be 


c 


0 


(2.31) 


where  g  is  gravitational  acceleration  and  k  is  the  wave  number.  The 
+  and  -  indicate  propagation  along  the  positive  and  negative  x  axes,  re¬ 
spectively  (the  x  axis  is  aligned  with  the  wave-number  vector  k) . 

In  the  presence  of  a  uniform  horizontal  current  whose  component  in 
the  direction  of  the  x  axis  is  UQ,  the  waves  are  Doppler  shifted  so 
that  the  phase  velocity  becomes 


c 


+  uo 


When  ocean  current  U  varies  linearly  with  depth,  c  becomes 


c  =  U 


z=-l/(2k) 


±  c0(l  +  s2) 


1/2 


where  s  =  U^/(2kc^).  Here,  denotes  the  current  component  in  the  x 

direction  and  IT  is  its  derivative  with  depth.  The  signs  are  associated 
with  the  propagation  direction  with  respect  to  the  x  axis. 

Comparison  of  the  second-order  solution  [Eq.  (2.30)]  to  the  above 
exact  expression  for  c  (for  linear  U)  suggests  a  scale  length  l* 
where 


l*  --  min 


41% 

U' 


X 


(2.32) 


for  checking  the  validity  of  the  perturbation  approach.  Here,  min  in¬ 
dicates  the  minimum  over  an  effective  depth  on  the  order  of  L.  To 
first  order  in  X.*/L, 
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} 


C 


±co  +  2k 


r° 

j  Ux(z)  exp(2kz)  dz 


(2.33) 


which,  from  the  definition  of  c,  can  be  expressed  in  terms  of  the  tem¬ 
poral  frequency  of  oscillation  f  =  0)/(2tt).  As  a  result, 


f  =  ±fB  + 


iMz)  exp(2kz)  dz 


(2.34) 


where 


f 


B 


!!!o  .iasi. 

2TT  2lT 


1/2 


is  the  frequency  of  oscillation  of  the  ocean  surface  wave  in  the  absence 
of  current. 

Roughly  speaking,  the  effective  Doppler  shift  is  caused  by  current 
in  the  uppermost  ocean  layer  having  a  thickness  on  the  order  of  one  wave¬ 
length,  which  is  evident  in  the  exponential  weighting  factor  exp(2kz). 
The  requirements  of  incompressibility  and  continuity  demand  a  divergence- 
free  wave-induced  velocity  field  for  the  fluid.  Because  this  velocity 
field  has  a  horizontal  scale  of  variation  equal  to  the  wavelength  L,  its 
vertical  scale  of  variation  must  also  be  equal  to  L.  This  being  true, 
it  is  influenced  by  the  ocean  current  within  a  depth  cn  the  order  of  one 
wavelength. 

The  ocean-current  magnitudes  encountered  in  this  investigation  are, 
at  most,  approximately  40  cm/sec,  and  still-water  wave  phase  velocity  is 
at  least  280  m/sec.  This  implies  that  the  effective  scale  length  defined 
in  Eq.  (2.31)  is  much  larger  than  the  vertical  scale  of  variation  of  U(z). 
As  a  result,  the  effects  to  second  order  in  i*/L  car.  be  ignored  in  most 
cases. 


23 


i 

* 

i 

? 

> 

I 


i 


Chapter  III 


STATISTICAL  DESCRIPTION  OF  OCEAN  WAVES 


The  purpose  of  this  chapter  is  threefold:  first,  to  discuss  the 
statistical  nature  of  the  ocean  waves  and  the  associated  statistical 
estimation  problem,  second,  to  introduce  the  concept  of  directional 
spectra  and,  third,  to  compare  their  various  definitions.  The  first 
two  are  of  direct  relevance  to  this  work;  the  third,  however,  requires 
clarification  and  should  be  useful  as  a  reference  for  radio  oceanogra¬ 
phers  . 

The  highly  irregular  nature  of  the  ocean  surface  makes  detailed 
measurements  difficult  to  obtain.  The  standard  approach  is  to  model 
the  irregular  temporal  and  spatial  variations  of  any  relevant  quantity 
pertaining  to  surface  motion  as  realizations  of  an  ergodic  random  pro-  j 

cess.  The  ensemble  average  of  a  function  of  this  quantity  can  be  as-  j 

! 

ymptotically  approached  by  averaging  the  function  over  an  increasing  j 

number  of  time  and  space  samples.  j 

* 

The  quantity  of  special  interest  is  the  long-crested  wave  compo-  ] 

nent  A(k,aj)  of  ocean-surface  motion  selected  by  the  Bragg-scattering  j 

mechanism.  If  N  is  the  elevation  of  the  surface  above  the  mean  sea 
level  z  =  0,  then  i 

A(k,u>)  =  jjj  dr  dt  n(r,t)  exp[i(k  •  r  -  wt)]  (3.1)  \ 

-00  I 

I 

where  | 

\ 

k  =  wave-number  vector  | 


a)  =  wave  radian  frequency 

r_  =  horizontal  spatial  vector  defining  the  spatial  location  of  any 
surface  element 

t  =  time 


This  is  a  divergent  integral,  and  the  problem  can  be  resolved  by  using 
the  Dirac  delta  function.  Note  that,  in  practice,  A(k,w)  is  a  quan¬ 
tity  that  cannot  be  determined  exactly  because  measurements  of  T)  for 


it' 


points  are  not  available;  at  best,  only  a  close 


all  temporal  and  spc.  ^  points  are  not  available;  at  best,  only  a  close 
estimate  of  A(k,w)  can  be  obtained. 

Section  A  reviews  the  problems  of  spectral  estimation  and  the  rela¬ 
tionship  between  ensemble  and  actual  averaging.  The  Wiener-Khinchin  the¬ 
orem  is  derived  in  Section  B,  and  various  power  spectra  used  by  oceanog¬ 
raphers  and  radio  oceanographers  are  compared. 


'll 

m 


A.  Practical  Spectral  Estimation 

According  to  Eq.  (3.1),  it  is  not  possible  to  measure  A(k,u>) .  This 
long-crested  wave  component  can  be  estimated,  however,  by  observing  the 
surface  elevation  p  over  finite  time  duration  T  and  finite  spatial 
extent  L. 

Without  neglecting  any  salient  characteristics  of  the  problem,  this 
discussion  is  limited  to  a  time  series  s(t)  whose  finite-time  transform 
at  time  X  and  frequency  f  is 


Sff 


,  i  rr/2 

=  5?  lT/2 


s(t  +  T)  w(t)  exp(-i2irft)  dt 


(3.2) 


where  w(t)  is  a  window  function  that  can  be  chosen  to  best  fit  this 
finite-time  transform  to  the  theoretical  transform 


S(f) 


.J-f 

2m  J 

«/_C 


s  (t)  exp(-i2lTft)  dt 


(3.3) 


Harris  [1978]  has  reviewed  the  many  chorees  of  w(t) .  In  the  discussion 
here,  a  rectangular  window  has  been  selected  with  width  T,  height 
unity,  and  centered  at  t  =  0,  and  Eq.  (3.2)  reduces  to 


S(f 


■T>  -  h  £ 


T/2 

T/2 


s (t  +  T)  exp(-i2mft)  dt 


(3.4) 


Because  the  phase  of  the  signal  s(t)  is  not  known  a  priori,  the 
above  Fourier-transform  kernel  does  not  adjust  itself  as  different  seg¬ 
ments  of  the  signal  are  examined.  There  are  the  following  two  major  dif¬ 
ferences  between  the  theoretical  transform  (3.3)  and  its  estimate  (3.4). 
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The  rectangular  time  window  leads  to  a  transform  that  is 
the  convolution  of  the  theoretical  transform  and  sin  (nfT)/ 
(fT)  whose  first  zeros  are  separated  by  2/T. 


•  The  estimate  can  be  time-dependent. 

A  A 

To  obtain  average  values  of  any  function  F(S)  of  the  estimate  S, 
S  and  F  are  computed  for  various  values  of  T  (t  =  nT  ,  n  =  0,1, . . . , 

^  cl 

N-l) ,  and  the  average  is  then  taken.  A  function  of  S  most  commonly 
used  is  its  squared  magnitude  (power  spectrum)  and,  to  simplify  the  fol¬ 
lowing  derivations,  the  function  F  is  assumed  to  be  S  itself.  Writ¬ 
ing  the  average  value  as  5(f), 


S(f) 


s  (t  +  nT  )  exp(-i2Tfft)  dt 
a 


(3.5) 


Bringing  the  summation  sign  through  the  integral,  it  is  apparent  that 

averaging  over  S  is  equivalent  to  averaging  over  the  time  signal 

s (t  +  nT  ) , 
a 

1  fT/2  - 

S (f )  =  —  I  s(t)  exp(-r2TTft'i  dt  (3.6) 

•'-T/2 


where 


s(t) 


N-l 

77  /  s(t  +  nT  ) 
N  4-  a 

n=0 


It  can  be  shown  that  the  averaging  process  is  equivalent  to  passing  s(t) 
through  a  filter  with  a  transfer  function  H(f)  such  that 

2 

sm  (irfNT  ) 

J- _ 5 

2  2 

N  sin  (TTfT  ) 
d 


H (f )  | 


and  this  filter  function  is  illustrated  in  Fig.  3.1.  If  the  average  sam¬ 
pling  time  T  is  small  so  that  the  sampling  rate  is  higher  than  twice 

ct 

the  bandwidth  of  the  signal,  this  discrete  averaging  process  is  equiva¬ 
lent  to  passing  the  signal  through  a  lowpass  filter  of  width  1/NT  . 

d 


Fig.  3.1.  MAGNITUDE  SQUARED  OF  FILTER  FUNCTION  VS 
NORMALIZED  FREQUENCY  f T  . 


Statistically,  the  average  over  discrete  samples  is  a  maximum-like¬ 
lihood  estimator  of  the  ensemble  mean  if  these  samples  are  obtained  from 
a  population  with  a  gaussian  distribution.  If  N  is  large  enough  for 
the  confidence  level  to  exceed  the  requirements,  this  estimator  is  con¬ 
sidered  the  ensemble  mean  and  is  denoted  by  angular  brackets  (  ).  Ap¬ 
plying  this  operator  to  Eq.  (3.4), 


exp  (~i2irft)  dt 


(3.7) 


because  s(t)  is  the  random  process.  Note  the  similarity  between  the 
time-average  operator  (denoted  by  the  overbar)  in  Eq.  (3.6)  and  the  en¬ 
semble-average  operator  in  Eq.  (3.7).  Both  operate  solely  on  the  random 
process  s(t). 

The  above  discussions  have  considered  the  problem  of  estimating  the 
ensemble  mean  in  practice.  In  the  following  analysis,  the  statistical 
approach  is  adopted  for  convenience. 


B.  Relationships  among  Various  Power  Spectra  Used  by  Radio 
Oceanographers 

In  this  section,  the  Wiener-Khinchin  theorem  is  derived  in  a  heuris¬ 


tic  manner,  based  on  the  following  assumptions. 


ws  ** ***fi»r<r*frhf*rir*r*** 


•  The  random  process  n(£»t)  has  an  ensemble  mean  h  that 
is  time-independent  and  spatially  uniform.  For  conveni¬ 
ence,  the  coordinate  system  is  adjusted  so  that  this  mean 
value  is  zero. 

•  The  autocovariance  function  defined  by 


C 

g 


V  11 ' -2 ' 


is  invariant  under  time  and  space  translation;  that  is. 


C  (r.  ,t^  ,-r^,t_)  =  C(r.  -  r„,t  -  t  ) 
g  — 1  1  —2  2  —1—2  1  2 


These  assumptions  in  the  time  (space)  domain  restrict  the  following  con¬ 
siderations  to  processes  that  are  stationary  (homogeneous)  at  least  to 
second  order,  and  a  power  spectrum  proportional  to  the  mean-squared  spec¬ 
tra]  amplitude  will  be  related  to  the  autocovariance  function  of  this 
process. 

Note  that  the  first  assumption  implies  that  the  spectral  amplitude 
defined  by  A(k,u))  is  also  zero  mean, 

(A(k,u»)  =  0 

The  covariance  of  spectral  amplitudes  at  points  (k,w)  and  (k* ,0)’)  is 
thus 


<A(k,0»  A*  (k '  ,(*}')) 

drdtdr’dt'  C(r-r',t-t')  exp[i  (k  •  r  -  k'  •  r  ’ )  -  i  (cot -U)’t' )] 

where  the  single  integral  represents  the  sixfold  integral  and  the  *  indi¬ 
cates  a  complex  conjugate.  The  variable  set  (r,r/ ,t,t')  is  now  trans¬ 
formed  to  (r,p,t,T)  where  p  =  r-r'  and  T  =  t-t’.  Recall  that 
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dt  exp(±iwt)  =  27r6(w) 


where  S(u)  is  the  Dirac  delta  function.  This  integral  is  only  conver¬ 
gent  in  the  sense  of  a  Cesaro  limit,  and  the  delta  function  has  a  symbolic 
meaning,  as  discussed  by  Born  and  Wolf  [1975].  Based  on  this  relation¬ 
ship,  the  covariance  function  above  reduces  to 

<A(k,0))  A*(k',0J’)) 


3  r00 

=  J  d£  dT  C'p,T)  exp[i(k  •  p  -  ojt)]  6(k  -  k ' )  6  (u>  -  to ' ) 


It  has  been  demonstrated,  therefore,  that  spectral  amplitudes  at 
different  wave  numbers  or  frequencies  are  uncorrelated.  This  result  as 
applied  to  a  wide-sense  stationary  time  series  has  been  discussed  in  de¬ 
tail  by  Helstrom  [I960].  A  power  spectrum  P(k,w)  is  now  defined  by 


P(k,u)  6(k  -  k')  6(w  -  0)’)  =  <A(k,w)  A* (k ',«’)) 


(3.8) 


from  which 


3  r°° 

P(k,w)  =  f~^)  I  d£  dt  C (£, T )  exp [i  (k  •  £  -  wt)] 

*"  '  •'—CO 


(3.9) 


is  obtained.  This  is  the  Wiener-Khinchin  theorem  wherein  the  power  spec¬ 
trum  defined  by  Eq.  (3.8)  is  the  Fourier  transform  of  the  autocovariance 

function  of  the  homogeneous  and  stationary  process.  From  the  inverse  re- 

2 

lation,  the  mean-squared  height  H  of  the  ocean  surface  is  determined 
to  be 


foo 

dk 

00 


du)  P (k,u>) 


(3.10) 


w'’-i 


From  the  definition  in  Eq.  (3.8),  the  power  spectrum  is  positive  real. 

In  addition,  because  the  autocovariance  function  is  real,  it  follows  that 


P(k,u>)  =  P(-k,-u>) 


(3.11) 
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For  first-order  gravity  waves  in  deep  still  water,  the  dispersion 
relation  has  only  two  roots  [Eq.  (2.21a)], 

0)  =  ±0)  =  ±  ✓gk" 

BO  a 


so  that  the  frequency  dependence  in  the  power  spectrum  P  can  be  removed, 
and  a  wave-height  directional  spectrum  S(k)  can  be  generated.  Because 
there  is  no  unique  definition  of  S  (k ) ,  however,  various  definitions 
found  in  the  literature  are  presented  in  sufficient  detail  for  direct 
comparisons. 

As  a  result  of  the  pioneering  work  of  Rice  [1951]  on  scattering  from 
random  surfaces ,  a  direct  extension  of  his  original  wave-number  spectrum 
W(k)  to  include  an  added  time  dimension  is  frequently  encountered.  Be¬ 
cause  his  spectrum  is  a  Fourier  analysis  in  space  only,  the  complex-con¬ 
jugate  relation  (3.11)  reduces  to 


W (k)  =  W(-k) 


(3.12) 


Generally,  P(k,oj)  is  not  symmetric  about  either  of  the  two  wave-number 

A  /\ 

axes  k^  and  k^.  To  generate  a  directional  spectrum  that  obeys  Eq. 
(3.12),  it  is  necessary  to  add  P  to  its  mirror  image  about  the  origin 
of  the  wave-number  axes.  The  new  directional  spectrum  (k) ,  there¬ 

fore,  is 


P(-k,w)  +  P(k,co)  =  ~  Srl(k)  ^6(0)  -  wB0)  +  6 


(w  +  (OH  (3.13) 


The  factor  preceding  s  ^ (k)  is  required  so  as  to  conform  with  the  nor¬ 
malization  relation  that  Rice  adopted, 


Joo 

dk 

00 


s  .  (k) 
rl  — 


(3.14) 


which  can  be  verified  by  substituting  Eq.  (3.13)  into  (3.10). 

To  depict  the  various  directional  spectra  simply,  the  following  ideal 
sea  conditions  are  assumed.  A  group  of  wave  trains  with  a  narrow  spread 
in  wavelength  is  propagating  in  the  +x  direction  while  another  group  with 
a  lower  energy  density  and  a  smaller  spread  in  wavelength  is  propagating 
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in  the  -x  direction.  The  exact  dependence  on  wave  number  k  is  assumed 
to  be  triangular — an  abrupt  cutoff  at  the  lower  wave-number  limit,  and  a 
linear  rolloff  from  there  to  the  higher  wave-number  limit.  The  directional 
spectrum  S^ (k)  along  the  k^  axis,  illustrated  in  Fig.  3.2,  describes 
the  ideal  sea. 


Fig.  3.2.  DIRECTIONAL  WAVE-HEIGHT  SPECTRUM 
Srl(k)  FOR  IDEAL  SEA  CONDITIONS.  The 
directional  spectrum  Sri (k)  is  shown 
only  in  the  (k^,0)  plane. 

The  inverse  relation  can  be  obtained  from  Eq.  (3.13)  as 

P(k,U)B0)  Aw  +  P(k,-wB0)  Aw  =  Srl(k)  (3.15) 

by  integrating  over  a  narrow  region  in  the  frequency  domain  centered  at 
w  =  wbq.  The  choice  of  Aw  (the  infinitesimal  frequency  interval  over 
which  P  is  constant)  is  determined  by  the  limiting  function  from  which 
the  delta  function  is  derived.  In  comparison  to  a  Fourier  series  in  the 
limit  of  infinite  period  T,  Aw  should  be  2tt/T.  From  Eqs.  (3.11)  and 
(3.15), 

S  _  (k)  =  S  (-k)  (3.16) 

rl  —  rl  — 

which  is  in  accordance  with  Rice's  constraint  [Eq.  (3.12)]. 

Because  there  is  an  awkward  factor  of  4  in  the  normalizing  relation 
in  Eq.  (3.14),  attempts  have  been  made  to  eliminate  it  by  introducing  a 
modified  spectrum  Sr2 (k)  as 

sr2(y  =  J  sri<3i>  (3.17) 
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From  this  definition,  the  volume  under  the  S  ^ (k)  surface  is  exactly 
equal  to  the  ocean  mean-squared  height. 

The  drawback  of  the  above  two  definitions  [Eqs.  (3.15)  and  (3.17)] 
is  that  the  energy  in  the  waves  propagating  in  the  +k  direction  is 
combined  with  the  energy  in  the  waves  propagating  in  the  -k  direction 
to  determine  the  height  of  S(k).  As  a  result,  S(k)  cannot  be  inter¬ 
preted  as  energy  associated  exclusively  with  waves  moving  in  the  4k  or 
-k^  direction.  To  eliminate  this  plus-minus  ambiguity,  Barrick  [1970] 
introduced  the  spectra  S+(k)  and  S  (k)  and  defined  the  x  axis  as 
the  direction  of  propagation  of  the  radio  wave  incident  on  the  ocean 
surface;  therefore, 

S±(k)  =  P(k,±oO  Au)  (3.18) 

where 

4o>bq  k  in  the  plus-x  half-plane 

-0)  k  in  the  minus-x  half-plane 

BO 


These  spectra  are  illustrated  in  Fig.  3.3  for  the  same  ideal  sea  condi¬ 
tions  described  for  Fig.  3.2.  This  is  an  extremely  awkward  definition 
because  the  energy  associated  with  waves  propagating  in  the  positive  k^ 
half-plane  is  contained  in  S+(k)  and  the  energy  for  negative  k^  is 


a.  S  (k)  in  the  (k  ,0)  plane 
+  —  x 


b.  S  (k)  in  the  (k^O)  plane 

Fig.  3.3.  DIRECTIONAL  WAVE-HEIGHT  SPECTRA 
FOR  IDEAL  SEA  CONDITIONS. 
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contained  in  S  (k) •  lc  is  not  clear  if  a  long-crested  wave  in  the  ±y 

direction  (for  which  k  =0)  should  be  included  in  S  or  S  . 

x  + 

A  simple  modification  of  Eq.  (3.17)  produces  a  directional  spectrum 
Sq  free  of  the  above  ambiguities, 

P(k,w)  =  r-  S  (k)  6(u  -  a>  )  +  q-  S  (-k)  5(u>  +  w  )  (3.19a) 

—  i.  o  —  B0  2o—  BO 

or,  equivalently, 

P(k,w  )  Ao>  =  q-  S  (k)  (3.19b) 

—  BO  2  o  — 

Again  for  the  same  ideal  sea  conditions  described  for  Fig.  3.2,  this  def¬ 
inition  of  the  directional  spectrum  is  illustrated  in  Fig.  3.4.  It  can 


Fig.  3.4.  DIRECTIONAL  WAVE-HEIGHT  SPEC¬ 
TRUM  SQ(k)  FOR  IDEAL  SEA  CONDITIONS. 
SQ(k)  is  shown  only  in  the  (kx,0)  plane. 


be  verified  that  the  complex-conjugate  relation  [Eq.  (3.11)]  is  not  vio¬ 
lated;  it  also  follows  that  the  mean-squared  height  is 

Joo 

dk  S  (k)  (3.20) 

00  “  °  “ 

Recalling  that  the  Fourier-transform  exponent  is  i  (k_  •  r_  •  ut) ,  Eq. 

(3.19b)  implies  that  SQ(k)  is  the  energy  density  associated  with  waves 
propagating  in  the  k_  direction  and  that  SQ(-k)  is  associated  with  the 
-k  direction.  The  definition  in  Eq.  (3.19a)  is  the  corrected  version 
of  that  obtained  by  Phillips  [1966,1978],  and  it  is  identical  to  that  of 
Longuet-Higgins  et  al  [1961]  except  that  the  inverse  of  their  power  spec¬ 
trum  is  the  complex  autocovariance  function  of  a  process  that  is  the  ana¬ 
lytic-signal  representation  of  the  surface-height  process  r)- 
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Fig.  3.5.  COMPARISONS  OF  VARIOUS  WAVE- 


HEIGHT  DIRECTIONAL  SPECTRA.  Numbers 
beside  the  shaded  area  indicate  heights 
of  pedestal  above  the  k^  -  k^  plane . 


Chapter  IV 


THEORETICAL  STUDIES  OF  FIRST-ORDER  BACKSCATTER 
FROM  OCEAN  SURFACE  WAVES  OF  DECAMETER  WAVELENGTH 


The  physical  processes  occurring  in  the  ocean  as  influenced  by  the 
presence  of  current  and  the  statistical  nature  of  the  ocean  surface  waves 
have  been  analyzed  in  Chapters  II  and  III.  In  this  chapter,  the  problem 
of  remotely  sensing  these  ocean  waves  by  means  of  an  HF  pulsed  radar  is 
discussed.  This  radar  system  consists  of  a  transmitter  and  a  receiver 
located  in  proximity  to  one  another  so  that  they  are  colocated  for  the 
theoretical  considerations  to  be  presented  here. 

As  shown  in  Fig.  4.1,  the  radar  system  is  situated  on  the  coast  at 

an  elevation  h  much  smaller  than  the  distance  rQ  from  the  ocean  patch 

it  is  designed  to  probe.  (In  the  experiments,  h  is  approximately  40  m 

and  rQ  is  roughly  10  to  40  km.)  To  a  good  approximation,  the  radio  wave 

transmitted  by  the  radar  propagates  at  a  zero  elevation  angle  with  respect 

to  the  sea  surface.  The  azimuthal  extent  L  of  the  ocean  patch  is 

Y 

established  by  the  beamwidth  of  the  radar  system  in  the  horizontal  plane; 
this  beamwidth,  in  turn,  is  defined  by  the  combined  beam  formed  by  the 


Fig.  4.1. 


RADAR  ILLUMINATION  OF  THE  OCEAN  SURFACE. 
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transmit  and  receive  antennas.  To  obtain  spatia]  resolution  in  the  radial 
direction,  the  radar  is  pulsed,  and  the  width  of  each  pulse  determines 
the  size  of  the  ocean  patch  in  the  radial  direction. 

Resonant  radio-wave  backscatter  from  the  ocean  surface  occurs  when 
the  surface  has  a  specific  spatial  structure;  its  wavelength  must  be  one- 
half  the  wavelength  of  the  radio  wave,  and  its  crests  must  be  transverse 
to  the  radio-wave  propagation  direction.  This  Bragg-scattering  phenome¬ 
non  has  been  explained  by  detailed  formulation  of  the  scattering  problem 
based  on  Rice's  perturbation  analysis  of  the  boundary  conditions  [Rice, 
1951]  and  the  Stratton-Chu  vector  diffraction  integral  [Barrick,  1972] . 

To  detect  ocean  current,  the  frequency  of  the  backscattered  signal  but 
not  its  absolute  amplitude  is  important.  The  following  derivations  based 
on  simple  wave-propagation  concepts  will  demonstrate  how  this  frequency 
is  related  to  ocean  current. 


A.  Physical  Theory  of  First-Order  HF  Backscatter 

Each  pulse  transmitted  by  the  radar  is  represented  by 

s  (t)  =  g(t  -  t  )  exp(i2irf  t) 
n  n  c 

where  n  signifies  that  the  pulse  is  initiated  at  time  t  =  tn'  and  g(t) 
is  the  pulse  envelope  with  the  rising  edge  at  t  =  0.  The  complex  expo¬ 
nential  is  the  analytic-signal  representation  of  a  sinusoidal  waveform  at 
carrier  frequency  f__,  and  its  phase  is  set  arbitrarily  to  zero.  Typi¬ 
cally,  many  cycles  of  the  sinusoid  are  contained  in  each  pulse  (a  few 
hundred  in  the  experiments) . 

The  pulsed  signal  will  propagate  along  the  radial  direction  £  in 
the  horizontal  plane  (see  Fig.  4.2)  and  will  scatter  in  all  directions 
when  it  encounters  the  undulating  sea  surface  denoted  by  z  =  h(r,t)  where 
z  is  the  vertical  ordinate  measured  upward  from  the  mean  sea  surface.  The 
strength  of  the  signal  component  scattered  back  toward  the  radar  must  be  a 
function  of  the  surface  slope  3t)/3r  (henceforth  denoted  by  n’)  along 
the  general  "look"  direction  of  the  radar.  For  small  T) ' ,  the  functional 
dependence  of  the  backscattered  signal  on  T]'  can  be  expressed  as  a  Tay¬ 
lor  series,  and  the  leading  term  in  this  series  expansion  is  TJ * .  To 
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first  order,  therefore,  the  backscat- 
tered  signal  strength  is  proportional 
to  n ' . 

In  Fig.  4.2,  the  angle  8  is 
measured  from  the  beam  axis  of  the 
radar,  which  is  also  the  x  axis. 

At  any  time  t,  the  signal  returning 
from  the  infinitesimal  ocean-surface 
patch  located  at  distance  r  and 
angle  8  is  proportional  to  the 
following  factors: 


RADAR  BEAM  CENTER, 
£ 


Fig.  4.2.  DEFINITION  OF 
THE  COORDINATE  SYSTEM. 


•  transmitted  signal  at  time  t-2r/c,  where  c  is  the 
speed  of  light  in  free  space,  which  takes  into  account 
the  round-trip  propagation  delay  of  the  radio  wave 

•  ocean-wave  slope  r} '  (:r,t  -  r/c)  ,  where  time  delay  r/c 
signifies  that  the  scattering  surface  is  located  dis¬ 
tance  r  away 

•  area  of  the  infinitesimal  ocean-surface  patch — namely, 
rdrdG 

•  combined  directivity  G(8)  of  the  transmit  and  receive 
antennas 

-2 

•  two-way  propagation  loss  r  experienced  by  the  radio 
wave 


The  received  signal  s<t)  is  the  sum  of  contributions  from  all  infini¬ 
tesimal  ocean-surface  patches  covered  by  the  radar. 


s(t) 


f17 

f"  / 

d8  G (8) 

1  rdr  n'(r,t 

'-TT 

J0  ‘ 

The  constant  of  proportionality  has  been  ignored. 

In  the  above  integral,  integration  over  r  actually  has  finite 

limits  as  the  result  of  the  finite  duration  of  the  pulsed  signal  sn(fc) 

whose  pulsewidth  is  denoted  by  T  .  Assume  that  s  (t)  is  zero 

pw  n 


this  same  signal  delayed  by  2r/c 


outside  the  time  interval  (0,T  ); 

pw 

will  then  have  the  following  property: 

exp  /i2irf  t  -  i2(Tf  — )  g  (t  -  t  -  — )  t  <  t  -  —  <  t  +  T 
\  c  cc/  \  nc/  n-  c  “  n  pw 

0  otherwise 

The  expression  for  the  received  signal  s  (t) ,  therefore,  simplifies  to 

s(t)  =  exp(i2ufct)  J  d6  G(6) 

exp(-i2te)  g(TDLy  -  (4.1) 

where 

Q 

rQ  =  —  (T  -  T  ^/2)  =  radial  range  of  center  of  interrogated 
p  ocean  patch  measured  from  the  radar 

Q 

L  =  —  T  =  radial  size  of  ocean  patch 
x  2  pw 

$  =  2irfc/c  =  free-space  radio-wave  wave  number 
=  t  -  t  =  fixed  time  delay 

As  is  evident  above,  the  fixed  time  delay  TDLY  •'■s  sampling  time  of 

s(t)  measured  from  time  t  at  which  the  most  recent  signal  s  (t)  was 

n  n 

transmitted.  The  ocean  patch  located  at  range  rQ  is  also  called  the 
radar  range  cell. 

Typically,  rQ  varies  from  tens  to  a  few  hundred  kilometers,  and 
the  linear  dimension  of  the  ocean  patch  is  from  a  few  to  a  few  tens  of 
kilometers;  therefore,  the  one-way  propagation  delay  r/c  is  on  the 
order  of  milliseconds  or  less.  Within  this  time  period,  no  significant 
variation  in  the  ocean-surface  slope  is  expected.  As  a  result,  T)’ (£/ 
t  -  r/c)  can  be  replaced  by  r)'(r,t),  and  the  slowly  varying  function 
r  ^  can  also  be  replaced  by  rQ^*  This  yields 
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r  +L  /2 

Ox  -1  .  _ 

r  dr  T] 

r-L  /2 
0  x 


c  rro+L 

s(t)  =  r  exp(i2irf  t)  I  d6  G(0)  I  X  dr  >)' (r,t)  exp(-i23r)  w(r) 

°  J  VV2 

(4.2) 

where  w(r)  =  g (TDLy  -  2r/c) . 

For  ease  of  interpretation  of  this  result,  the  circular  coordinates 
fr,0)  are  transformed  into  rectangular  coordinates  (x,y)  with  x 
pointing  in  the  beam-center  direction  (see  Fig.  4.2).  The  oeamwidth  is 
assumed  to  be  small  so  that 

r  =  x  (4.3) 

If  the  half-width  of  the  antenna  beam  is  0  the  azimuthal  extent  L 

0  y 

of  the  probed  ocean  patch  is  L  =  2r  sin  9„. 

y  0  0 

After  integration  by  parts  (see  Appendix  D) , 


r„+L  ./2 


s(t)  =  2ir.  exp(i2Trf  t) 
0  c 


I  0  x 

Jy  -T.  / 


dx  w(x) 


V2  w2 


G(y)  B  ’l(x,y,t)  exp(-i2gx) 


(4.4) 


If  antenna  directivity  G(y)  and  the  radial  weighting  function  w(x)  are 

assumed  to  be  constant,  this  twofold  integral  is  proportional  to  the  two- 

dimensionaJ  Fourier  transform  of  the  surface  height  at  the  wave-number 

set  (20,0)  convolved  with  a  two-dimensional  spectral  window  with  widths 

of  approximately  2TT/L  x  2ir/L  . 

x  y 

In  practice,  both  G(y)  and  w(x)  are  nonuniform  functions  and 
vary  relatively  slowly  with  their  respective  arguments.  Their  effects 
on  the  above  integral  are  the  same  as  the  effects  of  the  time  window  on 
the  Fourier  transform  of  a  time  series.  As  a  result,  resolutions  in  the 
spatial  frequency  domain  are  coarser  than  ^n/L^  x  27T/L  ;  the  degrada¬ 
tion  can  be  as  large  as  a  factor  of  2  [Harris,  1978,  for  example].  The 
radar  parameters  can  be  chosen  so  that  these  resolution  limits  are  much 
smaller  than  the  observable  spectral  width,  as  in  the  Stanford  four-fre¬ 
quency  radar  to  be  described  in  Chapter  V.  The  resolution  estimates 
2ir/Lx  and  2tt/L^,  however,  are  useful  order-of-magniiude  approximations. 
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In  the  limits  of  large  and  L  ,  the  integral  in  Eq.  (4.4)  rep¬ 

resenting  the  backscattering  process  becomes  an  infinitely  sharp  spatial 
filter  at  wave  number  (2(3,0) .  Only  one  specific  long-crested  ocean-sur¬ 
face  wave  component  is  selected  by  this  electromagnetic  process,  and  this 
component  has  a  wavelength  L  one-half  the  radio  wavelength  X  =  2tt/B 
and  propagates  along  the  radar  beam  axis.  Because  this  radio-wave  back- 
scattering  phenomenon  is  similar  to  Bragg  scattering  of  x-rays  from  crys¬ 
tals  [Brillouin,  1953] ,  it  is  also  called  Bragg  scattering. 

At  HF  (3  to  30  MHz) ,  the  Bragg-selected  ocean  waves  have  wavelengths 
in  the  decameter  range,  and  their  slope  magnitudes  are  smaller  than  unity 
even  under  fairly  rough  sea  conditions.  This,  in  turn,  implies  that  Bn 
is  small  and  that  the  first-order  effect  investigated  here  is  dominant. 
This  observation  served  as  the  basis  for  the  perturbation  calculations  of 
Rice  [1951] . 


In  the  expansion  of  r  in  Eq.  (4.3),  all  terms  of  order  higher  than 

y/r^  were  neglected.  The  most  important  effect  of  these  omitted  higher 

order  terms  is  in  the  phase  of  exp(-2iBr)  in  Eq.  (4.2).  An  estimate  of 

the  upper  bound  on  L  above  which  Eq.  (4.3)  is  not  valid  is  obtained  by 

limiting  the  maximum  phase  excursion  caused  by  the  second-order  term  to 

1/2 

tt/4 .  The  result  is  approximately  or.e-third  the  Fresnel-zone  size  (rQA) 

When  L  is  larger  than  the  Fresnel-zone  size,  the  above  theory  can 
be  modified  by  breaking  up  the  integration  in  0  in  Eq.  (4.2)  into  a  sum 

of  many  integrations,  each  over  a  small 

RADAR  range  of  0.  The  local  Cartesian  coordi- 

BEAM 

AXIS  nates  are  then  defined  for  each  of  these 


small  integration  zones  so  that  the  local 
x  axis  is  aligned  with  the  axis  of  the 
zone  (see  Fig.  4.3).  The  azimuthal  size 
of  each  zone  is  chosen  to  be  small  com¬ 
pared  to  the  Fresnel-zone  size  so  that 
the  above  theory  is  locally  valid  in  each 
integration.  The  modified  Bragg-scatter- 
ing  theory  then  states  that  the  backscat- 
ter  radar  is  sensitive  to  all  ocean  waves 
of  length  X/2  and  propagating  in  direc¬ 
tions  that  fall  within  the  antenna  beam. 


Fig.  4.3.  INTEGRATION  ZONES. 


When  Eq.  (4.4)  is  valid  and  the  width  of  the  backscattering  spatial 


filter  is  zero,  the  radar  will  detect  only  two  ocean  wave  trains — one 
propagating  toward  the  radar  along  its  beam  axis  and  the  other  propagat¬ 
ing  away  from  the  radar.  Both  wave  trains  have  a  wavelength  of  L=X/2, 
and  their  wave  height  can  be  represented  as 


T)(x,y,t)  =  B+(2g,0)  exp(i2iTf+t  +  i2(3x)  +  B  (2(3,0)  exp(i2irf  t  +  i2(3x) 

(4.5) 


Here,  the  subscripts  +  and  -  denote  wave  propagation  toward  and  away  from 
the  radar,  respectively,  and  B+  and  B  designate  the  respective  crest 
heights  of  the  wave  trains  above  the  mean  sea  level.  The  frequencies  f 
and  f  are  defined  in  Eq.  (2.34)  .is 


f±  »  «B  *  A 

„  UgS)172 

B  2ir 


A  = 


r° 

I  U  (z)  exp(4$z)  dz 

J-<x>  x 


(4.6a) 

(4.6b) 


(4.6c) 


As  discussed  in  Chapter  II,  f  is  the  frequency  of  oscillation  of  the 

B 

gravity  waves  in  the  absence  of  ocean  current,  \  is  induced  by  the  cur¬ 
rent  and,  because  the  wave  with  positive  frequency  is  propagating  in  the 
negative  x  direction  according  to  Eq.  (4.5) ,  is  the  current  compo¬ 

nent  in  this  same  direction  (pointing  toward  the  radar) . 

Substituting  the  expression  for  the  ocean  wave  height  in  Eq.  (4.5) 
into  (4.4)  and  ignoring  the  finite-resolution  effect  and  some  constants 
of  proportionality,  the  following  expression  for  the  received  signal  is 
obtained : 


s(t)  =  exp  i2ir(fc  +  A)  t 


][B+(2B,0)  exp(i2TTfBt)  +  B  (23,0)  exp (-i2irfBt) 


(4.7) 
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First-Order  Doppler  Spectrum 


The  Doppler  spectrum  of  the  received  signal  s (t)  is  determined  by 
taking  a  finite-time  Fourier  transform  [see  Eq.  (3.4)]  over  a  coherent 
integration  time  T  of  s (t)  and  forming  the  magnitude  squared  of  the 
result.  Here,  T  is  chosen  to  be  large  so  that  T  is  not  the  limit¬ 
ing  resolution  in  the  spectrum. 

In  the  absence  of  ocean  current.,  the  Doppler  spectrum  consists  of 

two  lines  (of  width  T  1)  located  at  f  +  f  and  f  -  f  (see  Fig. 

c  B  c  B 

4.4a).  Because  the  line  at  the  higher  frequency  is  the  result  of  Bragg 
scattering  from  ocean  wave  trains  propagating  toward  the  radar,  it  is 
referred  to  as  the  approaching  Bragg  line  and  f  as  Bragg  frequency. 
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Fig.  4.4.  THEORETICAL  DOPPLER  SPECTRA  OF  RADAR  RETURN. 


Based  on  the  same  reasoning,  the  line  at  the  lower  frequency  is  called 
the  receding  Bragg  line.  The  ratio  of  the  strength  of  the  first  Bragg 
line  to  that  01  the  second  is  the  ratio  of  the  mean-squared  wave  height 
of  the  approaching  wave  train  to  that  of  the  receding  wave  train. 

In  the  presence  of  ocean  current,  both  Bragg  lines  are  Doppler 
shifted  by  the  amount  A  which  is  proportional  to  the  ocean-current 
component  (in  the  direction  toward  the  radar)  averaged  over  depth 

with  an  exponential  weighting  function.  The  problem  of  detecting  ocean 
current  that  varies  with  depth  z  will  be  discussed  in  Chapter  VII.  To 
simplify  in  this  section,  is  assumed  to  be  uniform  with  depth,  and 

A  =  2U^/X.  The  displaced  Bragg  lines  are  illustrated  in  Fig.  4.4b. 

As  discussed  in  Section  A,  the  above  model  is  valid  if  L  is 

y 

smaller  than  the  Fresnel-zone  size.  It  has  been  demonstrated  that,  when 
L  is  many  times  this  limiting  size,  s(t)  is  the  sum  total  of  all 
contributions  of  backscatter  from  ocean  waves  of  wavelength  X/2  propa¬ 
gating  in  all  directions  that  fall  within  the  antenna  beam, 

s(t)  =  exp(i2irfct)  J  G(0)|b+  exp  |i2T:  (f^  +  A)  tj  +  B_  exp  j-i2ir  (fg  -  A)  t  j  d£ 

(4.8) 


Here,  A  =  -2r  •  U/X,  where  5?  is  the  unit  vector  along  the  radial  di¬ 
rection  from  the  radar.  The  current  vector  IJ  is  the  total  horizontal 
current  which  is  considered  uniform  with  depth.  The  ocean-wave  directional 

spectrum  (defined  in  Chapter  III)  is  assumed  to  be  omnidirectional  within 

2 

the  radar  beam  so  that  | B+ |  are  independent  of  0 . 

If  U  is  uniform  within  the  ocean  patch  probed  by  the  radar,  the 
shape  of  the  broadened  Bragg  lines  can  be  predicted  if  the  direction  of 
current  flow  is  known.  For  a  uniform  current  of  magnitude  UQ  flowing 
down  the  radar  beam  (in  the  negative  x  direction  in  Fig.  4.2),  the  cur- 
rent-induced  Doppler  frequency  A  will  take  on  values  ranging  from 
(2Uq/X)  cos  0q  to  2Uq/X,  where  0Q  is  the  half-width  of  the  radar 
beam.  Assuming  a  gaussian  beam,  the  Doppler  spectrum  is  illustrated  in 
Fig.  4.4c.  For  the  same  current  flowing  across  the  radar  beam  (in  the 
y  direction  in  Fig.  4.2),  A  then  takes  on  values  from  -(2Uq/X)  sin  0q 
to  (2Uq/X)  sin  0q.  The  resulting  spectrum  is  shown  in  Fig.  4.4d. 
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Observed  with  a  finite  radar  beamwidth,  the  maximum  width  of  each 

Bragg  line  attributed  to  U  is  (4UQ/A)  sin  9  .  If  the  ocean  current 

is  random  and  has  a  statistical  standard  deviation  0  ,  the  Bragg  width 

v 

resulting  from  this  random  component  is  on  the  order  of  20^/X. 

Two  other  aspects  of  the  ocean-surface  processes  that  can  cause 
broadening  of  the  Bragg  lines  are  the  coherent  length  of  the  long-crested 
waves  and  higher  order  ocean-wave  interactions. 


1 .  Effects  of  Finite  Coherent  Length  of  Long-Crested  Waves 


As  the  pulsed  radar  interrogates  the  surface  of  an  inhomogeneous 

sea  within  a  range  cell  of  finite  size  for  a  finite  duration  T,  long- 

crested  waves  with  different  heights  may  move  into  view.  For  waves  with 

1/2 

a  wave  number  of  23,  group  velocity  is  V  =  l/2(g/23) 

9 

required  for  them  to  move  through  the  entire  patch  of  length  L 


and  the  time 
is 


x 


t  = 
c 


L  /V  .  During  this  interval,  the  radar  essentially  samples  the  same  ocean 

patch  and  the  spectral  amplitude  B+  should  be  approximately  constant; 

that  is,  the  correlation  time  of  the  Bragg  amplitude  should  be  larger  than 

t  .  The  width  of  the  Bragg  line  is  therefore  narrower  than  t  1  if  t 
c  c  c 

is  smaller  than  the  coherent  integration  time  T.  Actually,  t"1  is  the 
limiting  resolution  of  the  finite-space  transform  because  is  the 

ratio  of  the  resolution  in  radian  frequency  to  the  resolution  in  wave 
number. 


Implicit  in  the  above  discussion  is  the  use  of  the  first-order 

dispersion  relationship  which  states  that,  for  a  particular  wavelength, 

there  can  be  only  two  frequencies  ±f  in  the  absence  of  current;  the 

B 

variation  of  the  Bragg-selected  wave  amplitude  with  time  is  determinis¬ 
tically  sinusoidal.  It  should  be  emphasized  that  this  dispersion  rela¬ 
tionship  is  valid  only  for  steady-state  conditions.  If  the  wind  changes 
significantly  within  the  coherent  integration  time  T,  an  additional 
temporal  variation  in  the  wave  amplitude  can  be  expected. 


2.  Higher  Order  Ocean-Wave  Interactions 

The  derivation  of  the  ocean-wave  dispersion  relationship  in 
Chapter  II  considered  only  the  terms  linear  in  the  wave  height  r).  Stokes 
(1847],  Longuet-Higgins  and  Phillips  [1962],  Johnstone  [j.975] ,  Huang  and 
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Tung  ]1976],  and  Weber  and  Barrick  [1977]  have  investigated  the  effects 
of  higher  order  terms  in  the  hydrodynamic  equations.  Some  of  their  rele¬ 
vant  results  are  summarized  in  this  section. 

Second-order  interactions  produce  a  dispersion  relationship 
0)  =  10  (k)  that  has  a  continuum  of  roots  at  the  specific  Bragg-selected 
wave-number  vector  3c  =  (23/0) .  Fortunately,  this  continuum  excludes  a 
region  centered  at  the  first-order  Bragg  frequencies  0)  =  ±/2g3  because 

of  the  low-frequency  cutoff  in  the  wave-height  spectrum  of  the  ocean.  The 
second-order  energies  thus  spread  themselves  in  a  continuous  frequency 
range  outside  the  first-order  Bragg  lines. 

Third-order  interactions  change  the  phase  speed  of  all  long- 

crested  waves,  including  the  Bragg-selected  wave  of  wave  number  (23,0). 

The  self-effect  (the  nonlinear  effect  in  the  presence  of  only  the  Bragg 

wave  itself)  is  an  increase  in  its  phase  speed  by  a  fraction  equal  to 
„  2 

1/2 (23a)  where  a  is  its  height.  This  is  called  "Stokes  current." 

Waves  with  various  wave  numbers  moving  in  parallel  or  opposite  to  this 
Bragg  wave  increase  or  reduce  its  phase  velocity,  respectively. 

Barrick  and  Weber  [1977c]  postulate  that,  because  this  phase- 
speed  change  is  dependent  on  the  ocean  wave-height  wh.ich  can  be  considered 
a  random  variable,  it  has  a  finite  variance  that  will  manifest  itself  as 
the  width  of  the  Bragg  line.  By  considering  only  the  effects  of  colinear 
waves  (long-crested  waves  moving  parallel  or  antiparallel  to  the  Bragg- 
selected  wave) ,  they  simplified  their  general  but  complicated  results  to 
a  tractable  form.  For  a  30  knot  wind  and  a  5  m  ocean  wave,  they  obtained 
a  phase-velocity  correction  of  5  percent  and  a  standard  deviation  (norma¬ 
lized  by  the  first-order  velocity)  of  0.05  percent.  This  correction  based 
on  mutual  interactions  will  have  a  sign  that  depends  on  the  wave  direction 
relative  to  the  wind;  waves  moving  with  the  wind  are  speeded  up,  and  waves 
moving  against  the  wind  are  slowed.  The  radar  cannot  distinguish  this 
correction  from  ocean  current.  The  0.05  percent  standard  deviation  is 
also  too  small  to  be  resolvable  by  the  radar. 


Chapter  V 


PESCADERO  EXPERIMENTS 


Since  May  1975,  ocean  backscatter  experiments  have  been  conducted 
at  a  site  located  near  Pescadero  on  the  northern  California  coast.  From 
May  through  December  1575,  the  experiments  were  run  approximately  four 
times  each  month,  and  a  typical  day  consisted  of  four  half-hour  runs, 
each  with  a  different  radar  pulse  width.  Similar  experiments  were  car¬ 
ried  out  during  June  1976  and  January  1977  with  nominal  radar  frequencies 
of  4.8,  6.6,  13.3,  and  21.7  MHz.  Based  on  the  decision  to  examine  sec¬ 
ond-order  interactions  at  higher  frequencies,  4.8  MHz  was  changed  to  29.8 
MHz.  Several  experiments  followed  during  August  and  September  1977  and 
January  1978  and,  at  that  time,  the  radar  pulse  width  was  fixed  at  50 
ysec  as  a  trade-off  between  the  signal-to-noise  ratio  and  spatial  reso¬ 
lution.  Approximately  four  hours  of  data  were  collected  daily,  and  con¬ 
current  in-situ  measurements  of  ocean-current  shear  were  obtained  by  de¬ 
ploying  spar  buoys  of  various  lengths  and  then  tracking  their  motions. 

The  Stanford  four-frequency  radar  system  used  for  data  collection 
is  described  in  Section  A.  The  data-processing  procedure  by  which  the 
sea-echo  Doppler  spectrum  was  obtained  and  the  subsequent  scheme  devel¬ 
oped  to  retain  only  the  first-order  Bragg  signal  are  discussed  in  Section 
B.  Possible  artifacts  introduced  by  the  radar  system  will  be  examined 
in  Section  C,  and  the  in-situ  measurements  are  analyzed  briefly  in  Sec¬ 
tion  D. 


A.  Stanford  Four-Frequency  Radar 

The  Stanford  radar  system  is  coherent  in  that  both  the  transmitter 
and  receiver  are  driven  by  the  same  free-running  stabilized  oscillator 
at  a  frequency  of  30  MHz.  Four  separate  phase-locked-loop  frequency  syn¬ 
thesizers  generate  the  transmitted  radar  frequencies  which  can  be  changed 
in  steps  of  10  kHz  from  the  front  panel  over  a  range  of  approximately  100 
kHz.  These  frequencies  can  be  adjusted  to  minimize  man-made  interference. 

The  transmit  antenna  is  a  vertical  half-rhombic,  approximately  250  m 
long  and  45  m  high  at  the  apex,  supported  by  a  helium-filled  balloon.  It 
is  located  on  a  flat  piece  of  land  40  m  above  sea  level.  The  vertical 
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plane  containing  the  half-rhombic  (referred  to  as  the  antenna  plane)  is 
in  the  direction  of  315 °T,  along  which  the  radar  system  is  about  1200  m 
from  shore.  The  antenna  wire  at  its  fixed  ends  has  roughly  1  m  of  slack, 
and  the  tension  is  taken  up  by  surgical  tubing;  motion  of  this  wire  in 
the  antenna  plane  is  thus  limited.  The  half-rhombic  is  supported  by  two 
guy  wires  transverse  to  the  antenna  plane.  Motion  in  this  transverse 
plane  is  again  limited;  however,  the  lack  of  elastic  support  permits 
larger  motions  here  than  in  the  antenna  plane.  The  advantage  of  this 
transmitter  system  is  its  relatively  high  gain  and  small  beamwidth  over 
a  wide  range  of  frequencies. 

The  receive  antenna  is  a  1  m  wideband  loop  located  50  yards  downhill 
from  the  transmit  antenna  in  the  direction  transverse  to  the  plane  of  the 
half -rhombic.  The  plane  of  the  loop  is  approximately  parallel  to  the 
plane  of  the  rhombic  so  that,  coupled  with  shielding  by  the  landscape, 
the  direct  effect  of  the  transmitted  pulses  is  minimized  although  not 
completely  eliminated  (see  Section  C) . 

In  Fig.  5.1,  the  measured  3  dB  beamwidth  of  the  combined  receive- 
transmit  system  is  plotted  vs  frequency.  At  5  MHz,  the  beamwidth  is 
approximately  25  degrees;  at  30  MHz,  it  is  roughly  12  degrees.  The  mea¬ 
sured  beam  center  moves  from  320°  T  at  30  MHz  to  290°  T  at  4  MHz  as  the 
result,  most  likely,  of  the  wavelength-dependent  guiding  effect  of  the 
surrounding  landscape  and  the  effective  attenuation  of  the  radio  wave  as 
it  propagates  over  land.  As  illustrated  in  Fig.  5.2,  the  radio  wave 
propagates  over  greater  stretches  of  land  on  the  north  side  of  the  an¬ 
tenna  beam  axis  than  it  does  on  the  south  side.  At  the  lower  frequency 
where  the  antenna  beamwidth  is  large  which  highly  accentuates  this  asym¬ 
metry,  the  direction  of  maximum  antenna  gain  thus  appears  to  have  moved 
south  from  315°  T. 


The  radar  is  pulsed;  a  trapezoidal  pulse  is  transmitted  every  400 
ysec.  The  nominal  pulse  width  T  (defined  as  the  width  between  the 
half-power  points)  is  selectable  from  the  front  panel  to  be  10,  25,  50, 
100,  or  200  ysec.  Four  different  carrier  frequencies  in  the  HF  band  of 
3  to  30  MHz  are  time  multiplexed  in  the  pulse  train  so  that  a  pulse  of 
the  same  frequency  repeats  every  1.6  msec  (see  Fig.  5.3).  The  peak  RF 
power  output  is  approximately  40  W. 
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The  receiver  system  is  shown  in  Fig.  5.4.  The  first  stage  of  ampli¬ 
fication  and  mixing  requires  four  channels — one  for  each  radar  frequency. 
The  local-oscillator  signal  for  each  channel  is  turned  off  during  a  blank¬ 
ing  period  lasting  through  the  width  of  the  transmitted  pulse  and  is  turned 
on  once  every  1.6  msec  when  the  preceding  carrier  frequency  is  30  MHz  below 
the  oscillator  frequency. 
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a.  Dual-conversion  superheterodyne  analog  system 
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b.  Digital  system 


Fig.  5.4.  RECEIVER  SYSTEM  FOR  THE  STANFORD  FOUR-FREQUENCY  RADAR. 
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The  signal  at  the  intermediate  frequency  is  then  passed  through  an 
amplifier  chain  and  a  crystal  filter  whose  width  is  approximately  40  kHz. 
The  second  conversion  stage  mixes  the  signal  to  the  final  intermediate 
frequency  of  0  MHz.  Two  parallel  mixers  driven  by  local  oscillators  that 
differ  in  phase  by  90°  are  used  to  distinguish  between  the  return  signal 
Doppler  shifted  below  and  above  the  transmitted  frequency.  The  two  chan¬ 
nels  are  designated  the  in-phase  channel  (I)  and  quadrature-phase  channel 
(Q). 

After  passing  through  identical  filters  and  amplifiers,  these  signals 
are  sampled  simultaneously  at  time  delays  of  typically  90,  130,  170,  and 
210  ysec  measured  from  the  leading  edge  of  the  transmitted  pulse.  Each 
delay  corresponds  to  t-t  ,  as  discussed  in  Chapter  IV.  Based  on  Eq. 
(4.1),  the  centers  of  the  four  range  cells  interrogated  are  at  9.75, 
15.75,  21.25,  and  27.75  km,  respectively,  from  the  radar,  assuming  a  pulse 
width  of  50  ysec. 

The  two  samples  at  each  range  bin  are  digitized  sequentially  by  an 
8-bit  A/D  converter.  Each  sample  is  then  stored  in  a  16-bit  accumulator 
designated  for  the  in-phase  or  quadrature -phase  channel  and  a  specific 
range  bin  and  radar  frequency.  There  are  two  sets  of  32  accumulators 
each;  in  each  accumulator,  N  (typically  140,  set  from  front  panel)  sam¬ 
ples  are  added  together.  While  one  set  is  waiting  for  N  time  samples 
to  accumulate,  the  contents  of  the  second  set  will  be  encoded  into  a  dig¬ 
ital  waveform  and  recorded  on  an  analog  tape  recorder.  Included  in  the 
data  stream  are  the  settings  of  all  the  front-panel  selector  switches  and 
an  internal  clock. 

The  summing  of  140  time  samples  lowers  the  sampling  rate  to  (140  x 

1.6  x  10“3)  ^  =  4.464  Hz  and  reduces  the  effect  of  white  noise  because 

the  operation  is  equivalent  to  a  lowpass  filtering  process.  The  filter 

,  -i 

function  is  periodic  with  period  (1.6  x  10  J)  =  625  Hz.  Each  period 
contains  a  sine  function  with  its  first  zeros  at  ±4.464  Hz,  and  the  sig¬ 
nal  falls  within  the  mainlobe  of  that  function.  The  effect  of  this  fil¬ 
ter  is  thus  to  reduce  the  white-noise  power  by  a  factor  of  140. 

If  the  model  in  Eq.  (4.7)  is  used  to  denote  the  signal  arriving  at 
the  receiver  front  end,  the  first-stage  and  final-stage  IF  signals  s^tt) 
and  s^  (t)  can  be  written  as 
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;^(t)  =  exp  J\ 2 rc ( 3  x  IQ7)  t  |b+  exp  Ji2Ti  (f^  +  A)  tj  +  B  exp  £-i.27T  (fg  -  A)  t  j 


+  complex  conjugate 


(5.1a) 


s2(t) 


=  B+  exp|i2T7(fB  +  A)  t  +  B  exp  — i27T  ( f ^  -  A)  t 


(5.1b) 


Here,  B+  and  B  are  the  amplitudes  of  the  ocean  wave  train  at  wave¬ 
length  A/2  propagating  toward  and  receding  from  the  radar,  respectively, 
and  A  is  the  current-induced  Doppler  shift  given  by  Eq.  (4.6c).  The 
in-phase  channel  is  assumed  to  be  the  real  part  of  a  complex  number  and 
the  quadrature  phase  channel  is  the  imaginary  part. 

From  Eq.  (4.6b) ,  the  Bragg  frequency  f  is 

£ 


f  =  0.102023  /f~  Hz 
B  c 


(5.2) 


where  f^  is  the  radar  carrier  frequency  in  megahertz.  The  values  of 
fB  for  the  five  nominal  radar  frequencies  used  in  the  experiments  are 
tabulated  in  Table  5.1,  together  with  the  wavelength  and  phase  velocity 
of  the  Bragg-resonant  ocean  surface  waves. 


Table  5.1 

CHARACTERISTICS  OF  BRAGG-SELECTED  OCEAN  WAVES. 

Bragg  frequency  f  is  from  Eq.  (5.2). 

B 


Radar 

Frequency 

(MHz) 

Ocean-Wave 

Wavelength 

(m) 

Phase 

Velocity 

(m/sec) 

Bragg 

Frequency 

(Hz) 

4.8 

31.25 

6.99 

0.224 

6.7 

22.39 

5.91 

0.264 

13.3 

11.28 

4.20 

0.372 

21.7 

6.91 

3.29 

0.475 

29.8 

5.03 

2.80 

0.557 

•&V.  '3 


B.  Data  Reduction 


To  obtain  the  first-order  signal  from  the  recorded  data,  these  data 
are  reduced  through  the  following  steps. 

Step  1 

The  recorded  signal  is  played  back  through  a  PCM  decoder  that  recov¬ 
ers  the  clock  signals  ana  converts  the  serial  data  stream  into  8-bit  data 
words.  These  data  words  include,  for  example,  front-panel  settings,  radar 
frequencies,  and  pulse  width.  The  decoded  data  are  stored  on  a  nine-track 
magnetic  tape. 

Step  2 

Depending  on  the  frequency  resolution  desired,  a  corresponding  num¬ 
ber  of  data  samples  are  extracted  from  the  magnetic  tape  for  each  range 
bin  and  carrier  frequency  and  are  stored  temporarily  as  16  sets  of  data 
on  the  disk.  The  observed  current-induced  Doppler  shift  ranges  from  a 
few  cm/sec  to  40  cm/sec.  To  resolve  the  small  current  and  still  retain 

reasonable  time  resolution  requires  a  coherent  integration  time  of  917.504 

12 

sec  which  corresponds  to  2  =  4096  complex  pairs  of  time  samples  for 

each  Fourier  transform.  Because  of  computer-memory  limitations,  the 
number  of  samples  is  reduced  by  a  factor  of  2  by  replacing  every  tv/o 
adjacent  time-sample  pairs  by  their  sum;  the  effective  sampling  rate 
therefore  becomes  2.232  Hz.  The  data  sets  are  then  retrieved  from  the 
disk  storage,  one  set  at  a  time,  and  the  following  processing  is  per¬ 
formed  on  the  selected  data  set. 

Step  3 

The  temporary  loss  of  synchronization  in  step  1  will  result  in  erro¬ 
neous  data  values.  The  frequency  of  occurrence  of  this  event  is  minimized 
by  selecting  a  unique  synchronization  pattern.  When  this  loss  occurs,  it 
affects  only  the  associated  data  frame  consisting  of  16  time-sample  pairs — 
one  for  each  range  bin  and  each  radar  carrier  frequency.  To  remove  such 
data  spikes,  an  algorithm  is  devised  whereby  the  mean  s  and  standard 
deviation  0  of  the  channel  (in-phase  or  quadrature-phase)  are  computed, 


*n  v-*-  <*v<i*»j.*«w'***j**  "- 


the  data  in  tnat  channel  are  searched,  and  all  data  points  that  deviate 
from  the  mean  s  by  more  than  3a  are  eliminated.  In  this  step,  the 
mean  is  removed  from  each  channel. 

Step  4 

A  Hamming  window  is  then  applied  to  the  in-phase  and  quadrature- 
phase  channels  in  the  data  set  of  2048  complex  pairs.  This  window  pro¬ 
duces  a  spectral  resolution  of  approximately  one-half  that  of  a  rectan¬ 
gular  window  (2.2  millihertz) ;  however,  the  spectral  sidelobes  are  lower 
than  that  of  the  rectangular  window  by  30  dB  [Harris,  1978] . 

Step  5 

This  set  of  windowed  data  is  then  Fourier  transformed,  using  the 
Cooley-Tukey  algorithm  of  the  fast  Fourier  transform  (see,  for  example, 
Rabiner  and  Gold,  p.  367) .  The  power  (Doppler)  spectrum  is  obtained  by 
taking  the  sum  of  the  squares  of  real  and  imaginary  parts  of  the  trans¬ 
formed  data.  The  results  are  stored  on  a  nine-track  magnetic  tape. 

Step  6 

The  processing  steps  3  through  5  are  repeated  for  each  of  the  data 
sets  stored  on  the  disk. 


Step  7 

Power  spectra  are  obtained  at  time  increments  of  one-half  the  coher¬ 
ent  integration  time  (approximately  7-1/2  min) .  This  is  possible  because 
of  the  50  percent  overlap  in  the  data  used  in  the  transform;  that  is,  the 
latter  half  of  the  data  used  in  obtaining  the  power  spectrum  of  a  previous 
time  interval  becomes  the  first  half  of  the  data  for  the  present  time  in¬ 
terval.  Preferential  weighting  of  the  Hamming  window  suppresses  informa¬ 
tion  contained  at  the  outer  edges  of  the  window,  and  the  50  percent  time 
overlap  recovers  this  lost  information. 
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Step  8 


The  large  variability  (or  jaggedness)  of  the  individual  power  spectra 

indicates  the  randomness  of  the  ocean-surface  process.  In  fact,  the  power 

spectral  estimate  behaves  as  a  chi-square  with  four  degrees  of  freedom 

(see  Appendix  A) .  To  reduce  this  variability,  a  smoothed  spectrum  is 

obtained  by  convolving  it  with  a  gaussian  window  whose  width  at  1/e  (or 

4.34  dB)  is  chosen  to  be  proportional  to  the  square  root  of  the  carrier 

frequency.  This  frequency  dependence  is  based  on  the  observed  behavior 

of  the  Bragg  width  (see  Chapter  VIII).  In  terms  of  frequency  bins,  each 

1.09  ird  llihertz  in  size,  this  window  width  is  8  bins  at  4.8,  10  bins  at 

6.8,  13  bins  at  13.3,  17  bins  at  21.8,  and  21  bins  at  29.75  MHz.  Because 

the  original  unsmoothed  spectrum  is  used  in  the  actual  determination  of 

the  Bragg-line  position,  these  smoothing  widths  are  not  critical. 

In  the  typical  smoothed  spectra  in  Fig.  5.5,  narrow  bands  of  strong 

signals  can  be  observed  at  frequencies  to  within  a  few  percent  of  ±f  . 

B 

As  discussed  in  Chapter  IV,  the  band  of  signals  near  +f  is  energy  scat- 

B 

tered  from  ocean  waves  propagating  toward  the  radar  and,  for  brevity,  are 

called  the  approaching  Bragg  line;  che  band  near  -f  is  called  the  re- 

B 

ceding  Bragg  line  for  similar  reasons.  In  the  immediate  neighborhood  of 
each  band,  a  continuum  of  lower  level  signals  attributed  to  second-order 
electromagnetic  and  hydrodynamic  effects  iTyler  et  al,  1974;  Barrick, 

1972;  Johnstone,  1975]  is  frequently  found  and,  most  often,  it  is  sepa¬ 
rated  from  the  first-order  band  by  a  sharp  null.  On  days  of  relatively 
calm  sea  and  particularly  at  the  farthest  range  cells,  this  second-order 
continuum  is  completely  submerged  in  system  noise  which  is  approximately 
white.  In  either  case,  further  analysis  of  the  first-order  bands  is  nec¬ 
essary,  and  this  is  made  possible  by  a  frequency-domain  filtering  scheme 
that  can  be  applied  to  either  first-order  band. 

Step  9 

Examination  of  the  Doppler  spectra  derived  from  the  entire  Pescadero 
data  set  reveals  that  the  current-induced  Doppler  shift  corresponds  to  an 
effective  uniform  current  having  a  magnitude  of  no  more  than  50  cm/sec. 
This  fact  is  used  in  the  search  for  the  peaks  of  the  Bragg  signals. 
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Fig.  5.5.  SMOOTHED  DOPPLER  SPECTRA  AT  FOUR  RADAR  FREQUENCIES  AND  FOUR 
RANGE  CELLS.  Data  were  collected  on  25  January  1978. 


The  purpose  of  the  frequency-domain  filtering  scheme  is  to  determine 
the  cutoff  points  on  either  side  of  the  band  peak.  An  algorithm  has  been 
developed  that  looks  at  the  local  slopes  of  the  logarithm  of  the  smoothed 
spectrum.  Each  is  computed  as  the  slope  of  the  linear  fit  to  a  number  of 
spectral  points,  and  the  number  is  the  same  as  that  for  the  smoothing 
window.  On  one  side  of  the  peak,  the  slope  must  be  positive  and  on  the 
other  side  it  should  be  negative.  This  scheme  establishes  the  point  on 
the  positive-slope  side  at  which  the  local  slope  becomes  smaller  than  a 
positive  preselected  cutoff  value  and  the  point  on  the  negative  -  slope 
side  where  the  local  slope  becomes  larger  than  a  negative  cutoff  value. 
These  two  points  then  define  the  limits  of  the  first-order  band.  The 
magnitude  of  the  cutoff  value  is  0.3  dB/millihertz. 
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Note  that  the  smoothed  spectrum  is  used  to  obtain  the  cutoff  points. 
These  can  then  be  considered  as  the  iuwci  and  upper  frequency  cutoffs  of 
a  rectangular  filter  to  be  applied  to  the  original  unsmoothed  spectrum; 
that  is,  they  define  the  first-order  band  in  the  unsmoothed  spectrum. 
Figure  5.6  is  a  block  diagram  of  this  process. 


FIRST- 

ORDER 

SIGNAL 


Fig.  5.6.  RETRIEVAL  OF  FIRST-ORDER  SIGNAL  FROM  DOPPLER 
SPECTRUM  OF  SEA  ECHO. 


C.  System-Introduced  Artifacts 

1.  Effects  of  Transmit-Antenna  Motion 

A  disadvantage  of  a  balloon-supported  antenna  is  that  its  mo¬ 
tions  may  broaden  the  spectral  lines.  To  understand  this,  the  following 
extension  of  Eq.  (5.1b)  is  used  as  a  model  for  the  postdetection  signal: 

s2(t)  =  Ap(t)js+  exp|i2it(fB  +  A)  tj  +  B_  exp |^-i2 n (f B  -  A)  t 
+  AQ(t)  +  AL(t)  +  D 
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(5.3) 


Hers,  A^(t)  is  the  transmitted  signal  whose  amplitude  slowly  varies 

with  time  because  of  the  motions  of  the  antenna  in  the  vertical  plane 

containing  the  beam  axis.  The  receive  antenna  is  inevitably  sensitive 

to  the  tail  of  the  transmitted  pulse  that  continues  through  the  receive 

period.  This  results  in  a  detected  direct  signal  Ap(t)  whose  time 

variation  is  caused  by  the  motion  of  the  half-rhombic  in  a  vertical 

plane  transverse  to  the  beam  axis.  This  signal  is  expected  to  decay 

with  range  as  the  transmit  pulse  tapers  off  with  time.  The  component 

A  (t)  is  caused  by  the  direct  reflection  from  stationary  targets;  its 
L 

temporal  variation  follows  the  antenna  motion  along  the  line-of -sight 

path  to  the  stationary  target.  The  constant  term  D  is  entirely  the 

result  of  leakage  of  the  30  MHz  local-oscillator  frequency  into  the 

signal  path  within  the  radar  box;  this  appears  as  a  line  in  the  Doppler 

spectrum  at  zero  frequency,  and  its  magnitude  is  several  orders  higher 

than  A  (t)  and  A  (t) . 

D  L 

The  general  modulation  function  A(t)  is  proportional  to 


}  exp [j2irL(t)/X] ,  where  L(t)  is  the  time-varying  component  of  the  path- 

«• 

*  length  from  the  source  to  the  receiver,  and  X  is  the  radio  wavelength. 

j  If  L (f.)  has  amplitude  L  and  frequency  of  oscillation  f  then,  for 

small  LQ/X,  the  spectrum  of  A(t)  consists  of  lines  with  amplitudes 

on  the  order  of  (L  /X)n  at  frequencies  of  f  =  ±nf  ,  n  =  0,1,...,°°. 
j  0  n  A 

|  For  actual  random  antenna  motions,  the  spectral  width  of  A(t)  is  roughly 

|  (L/X)  •  1/T  ,  where  L  is  the  average  amplitude  and  T  is  a  represen- 

|  tative  time  scale.  The  effect  of  antenna  motions,  therefore,  increases 

I  with  radar  frequency. 

|  The  resulting  Loppler  spectrum  from  Eq.  (5.3)  will  consist  of 

|  bands  of  energy  at  ±f  and  at  dc  where  width  is  dictated  by  transverse 

i  B 

‘  motions  with  higher  amplitude,  as  discussed  in  Section  A;  this  is  very 

\  apparent  at  higher  radar  frequencies  under  moderate -to-high  wind  condi¬ 


tions  (see  Fig.  5.5).  Note  that  the  effect  of  D  has  been  removed.  At 
±f  ,  there  is  a  convolution  of  the  antenna  motion-induced  spectrum  and 
the  ocean  wave-induced  and  current-induced  spectrum.  If  the  system  arti¬ 
fact  is  dominant,  the  bands  at  ±f  are  expected  to  exhibit  structures 

B 


•  The  structures  of  the  spectral  bands  at  ±fg  must  be 
similar. 

•  Band  structures  must  not  vary  with  range  from  the  radar. 


In  Fig.  5.7,  the  fine-scale  structures  of  the  approaching  and 
receding  Bragg  lines  are  shown  by  plotting  the  portions  of  the  unsmoothed 
spectra  centered  at  ±f  .  Note  that,  at  either  of  the  two  radar  frequen- 
cies,  the  spectra  are  dissimilar  at  different  ranges  and  that  the  ap¬ 
proaching  and  receding  lines  do  not  exhibit  the  same  characteristics.  For 
example,  at  21.75  MHz,  the  widths  of  the  receding  Bragg  lines  are  signif¬ 
icantly  larger  than  those  of  the  approaching  lines  and,  at  29.75  MHz,  the 
width  of  the  approaching  lines  at  the  near  range  is  much  smaller  than 
that  at  the  farther  range.  Such  dissimilarities  between  the  receding  and 
approaching  Bragg  lines  and  between  those  at  various  range  cells  are 
continuously  observed  in  the  Doppler  spectra  obtained  from  the  Pescadero 
data.  As  a  result,  it  is  unlikely  that  antenna  motions  in  the  plane  of 
the  half-rhombic  are  large  enough  to  increase  significantly  the  width  of 
the  Bragg  lines. 

Even  if  antenna  motions  are  large  enough  to  increase  the  width 
of  the  Bragg  line,  they  cannot  change  the  mean  position  of  that  line 
because  they  are  zero-mean  random  motions. 

2.  Local-Oscillator  Instability 


It  has  been  suggested  that  the  spectral  width  of  the  land  echo 
is  a  measure  of  the  local-oscillator  (LO)  stability.  This  echo,  reflected 
from  a  stationary  object,  is  at  the  frequency  of  the  transmitted  signal. 

It  thus  appears  at  zero  frequency  (dc)  after  IF  conversion,  and  its  width 
is  determined  by  the  stability  of  the  local  oscillator  over  its  time  of 
flight  which  is  on  the  order  of  tens  of  microseconds. 

As  evident  from  Eqs.  (5.1b)  and  (5.2),  the  Bragg  signal  will  be 
affected  by  the  variation  of  the  oscillator  frequency  over  the  integration 
time  of  15  min  required  to  produce  the  Doppler  spectra.  By  differentiating 
Eq.  (5.2),  it  can  be  shown  that  the  percentage  change  in  the  carrier  fre¬ 
quency  is  twice  that  in  the  Bragg  frequency.  Typical  values  for  the  Bragg 
frequency  and  Bragg  width  are  500  and  10  millihertz,  respectively.  If  this 
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4. 


Finite  Range-Cell  Size 


As  discussed  in  Chapter  IV,  the  finite  range-cell  size  x 

implies  that  the  width  of  the  two-dimensional  spatial  Bragg  filter  is 

Ak  x  Ak  where ,  to  a  good  approximation ,  Ak  =  2tt/L  and  Ak  =  2tt/L  . 
xy  xxyy 

To  translate  this  spatial-frequency  width  into  Af  =  Au)/2Tr  in  the  fre¬ 
quency  domain,  a  small  change  Ak  in  the  wave  number  is  related  to  a 
small  change  Aw  in  the  wave  frequency  through  the  wave  group  velocity 

V  , 

9 

Aw  =  (Ak)  V 

g 


where  Ak  is  the  magnitude  of  the  vector  with  components  Ak^,Ak  .  If 


/  -2  -2 \ 

L  =  (Lx  +  Ly  ) 

the  expression  for  Af  is 

V  .  V 

=  (5.5) 


where  is  the  wave  phase  velocity  tabulated  in  Table  5,1.  The  worst- 
case  Af  occurs  at  the  lowest  radar  frequency  of  4.8  MHz  and  the  short¬ 
est  pulse  width  of  25  ysec.  For  a  distance  of  15  km  and  a  beamwidth  of 

30°  (see  Fig.  5.1),  L  «=  7.5  km  and  L  3.8  km  so  that  f  «  1  mil- 

y  x 

lihertz.  This  worst-case  Af  is  slightly  smaller  than  the  resolution 
cell  of  the  918  sec  transform.  The  effect  of  finite  range-cell  size, 
therefore,  is  not  observable  in  the  data. 


D.  In-Situ  Measurements  of  Ocean  Current  by  Spar  Buoys 

Concurrent  in-situ  drifter  measurements  of  the  ocean  current  were 
obtained  during  several  of  the  radar  experiments.  The  drifters  were  four 
lengths  of  spar  buoys  built  of  plastic  pipes  loaded  with  an  appropriate 
amount  of  lead  shot  so  that  they  immersed  upright  in  water  to  predeter¬ 
mined  depths  of  1,  3,  6,  and  12  m,  respectively. 

The  locations  of  the  buoys  were  determined  by  a  microwave  tracking 
system  whose  range  resolution  was  on  the  order  of  a  few  meters.  As 
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illustrated  in  Fig.  5.8,  two  microwave  transponder  (or  range)  stations 
were  located  at  Pillar  Point  and  Pescadero,  respectively.  A  boat  used 
to  deploy  the  buoys  also  housed  the  master  station  that  continuously 
interrogated  the  two  transponders  on  the  coast  and  determined  the  dis¬ 
tance  of  the  boat  from  each  range  station  by  monitoring  the  round-trip 
time  delay  of  the  microwave  signals.  The  position  of  each  buoy  was  mea 
sured  by  carefully  pulling  the  boat  to  within  a  few  meters  of  it  and 
then  noting  the  distances  displayed  on  the  tracking  system.  These  mea¬ 
surements,  taken  at  time  intervals  of  an  hour  or  longer,  determined  the 
mean  drift  velocity  of  each  buoy.  Because  this  drift  was  typically  on 
the  order  of  1  km,  the  uncertainty  introduced  was  less  than  1  percent. 


Fig.  5.8.  LOCATIONS  OF  SPAR  BUOYS  AND  RANGE  STATIONS. 


The  buoys  were  not  completely  immersed  in  water  so  that  they  would 

be  visible  to  the  human  eye.  The  length  SL  remaining  above  the  sea 

surface  is  approximately  25  percent  or  the  total  length  for  the  shortest 

buoy  and  10  percent  for  the  longest.  The  effect  of  wind  drag  on  the  drift 

velocity  of  the  buoys  is  estimated  as  follows.  The  component  of  wind 

in  the  direction  of  current  flow  is  denoted  by  U  ,  current  speed  by  U  , 

a  w 

air  and  water  densities  by  p  and  p  ,  respectively,  and  the  immersed 

a  w 

buoy  length  by  Z  .  At  equilibrium,  U  <  U  <  U  and,  because  the  ac- 
w  w  D  a 

celerating  force  from  the  wind  is  balanced  by  the  retarding  force  from  the 
water , 
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The  ratio  p  /p  is  1.2  x  10  . 
w  a 

shortest  buoy  where  i  /£  =  1/3. 

w  a 

current,  therefore,  is 


The  effect  of  wind  drag  is  largest 

The  buoy  drift  U  relative  to 
Dr 
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where  is  assumed  to  be  small  compared  to  U^.  The  maximum  wind  speed 

at  a  10  m  height  was  approximately  5  m/sec  (10  knots)  during  the  January 
1978  experiments.  At  the  elevation  of  the  exposed  buoy,  the  wind  speed 
is  expected  to  be  lower  by  a  factor  of  2.  The  maximum  UDr  is  thus  -4 
cm/sec  which  is  smaller  than  the  statistical  uncertainty  in  the  buoy-in¬ 
ferred  drift  (see  Chapter  VII) . 

It  will  be  demonstrated  in  Appendix  B  that,  in  a  current  that  varies 
with  depth,  a  buoy  immersed  to  depth  l  will  drift  with  a  speed  equal  to 
the  current  at  depth  al,  where  a  =  0.5  and  0.27  for  current  varying 
linearly  and  logarithmically  with  depth,  respectively.  With  four  buoy 
lengths,  current  can  be  probed  at  four  different  depths  if  it  has  one  of 
the  above  profiles.  There  is  no  general  solution  to  the  nonlinear  prob¬ 
lem  of  determining  the  ocean-current  profile  from  the  buoy  drifts. 

To  obtain  some  indication  of  the  variation  in  ocean  current  over  a 
radar  range  cell,  two  sets  of  buoys  separated  by  approximately  5  km  were 
deployed  within  the  third  range  cell  (see  Fig.  5.8).  Because  only  one 


boat  was  used,  the  measurements  of  current  at  these  two  points  were  also 
separated  by  the  time  required  by  the  boat  to  move  from  one  location  to 
the  next  and  to  search  for  the  buoys — typically,  30  to  60  min. 


Chapter  VI 


DETERMINATION  OF  AVERAGED  OCEAN  CURRENT  FROM  DOPPLER  SPECTRA 


In  Chapter  V,  a  procedure  was  described  by  which  the  Doppler  spec¬ 
trum  of  the  sea-returned  signal  was  obtained.  An  algorithm  was  developed 
to  determine  the  spectral  points  associated  with  first-order  backscatter. 
These  spectral  points  from  the  first-order  band  at  either  the  positive  or 

negative  Bragg  frequency  f  are  used  in  this  chapter  to  determine  the 

B 

averaged  current  at  the  various  range  cells  probed  by  the  HF  radar  at  four 

frequencies.  The  results  obtained  in  the  previous  chapters  indicated  that 

the  first-order  band  at  the  positive  Bragg  frequency  is  Bragg-scattered 

energy  from  ocean  waves  moving  toward  the  ocean.  When  this  energy  band 

is  of  infinitesimal  width,  its  displacement  A  from  position  f  has 

B 

been  shown  [Eq.  (4.6b)]  to  be 


2tt 

26 


A 


46 


f 

•'n 


Vz) 


exp(-46z)  dz 


(6.1) 


where  U^(z)  is  the  component  of  current  in  the  direction  toward  the 
radar  and  z  is  depth  below  the  ocean  surface.  The  quantity  on  the 
left-hand  side  is  called  the  current-induced  Doppler  velocity  LJ  ( p) . 
Introducing  the  nondimensional  variable  z  =  46z,  Eq.  (6.1)  can  be  re¬ 
written  as 


U(6) 


exp (-z) 


dz 


(6.2) 


The  measured  Doppler  velocity,  therefore,  is  a  weighted  average  of  the 
actual  current  profile,  and  it  varies  with  the  probing  radio  wave  number 
because  the  weighting  function  has  a  scale  length  inversely  proportional 
to  6. 


As  has  been  demonstrated  in  Chapter  V,  the  observed  Bragg  line  is 
always  of  finite  width  which  is  assumed  to  be  a  manifestation  of  the 
turbulent  character  of  the  current  field;  that  is,  current  varies  nonde- 
terministically  both  temporally  and  spatially  within  the  radar  range 
cell.  This  assumption  is  experimentally  verified  in  Chapters  VII  and  VIII. 
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The  centroid  of  the  first-order  energy  band  will  be  used  as  an  estimate 

of  the  position  of  the  Bragg  line,  and  its  deviation  from  f  is  an  es- 

B 

A 

timate  of  A.  The  Doppler  velocity  U($)  thus  derived  will  be  shown  to 
be  the  space-time  average  of  the  underlying  current  profile. 

A  similar  argument  applies  to  the  first-order  band  at  -f  .  A  space- 
time  average  of  the  current  profile  can  also  ba  obtained  from  energy  scat¬ 
tered  by  ocean  waves  receding  from  the  radar. 


A.  Estimator  of  Current-Induced  Doppler  Velocity 


The  centroid  of  the  first-order  band  is  applied  in  this  section  to 

estimate  the  position  of  the  Bragg  line.  Specifically,  if  m  and  m 

U  J-» 

are  the  upper  and  lower  cutoff  points  as  determined  from  the  algorithm 
in  Chapter  V,  the  Bragg  position  is 


f 
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(6.3) 


where  the  sums  are  taken  from  m  =  m^  to  m  =  my,  and  is  the  spec¬ 

tral  value  at  the  Doppler  frequency  f  .  The  current-induced  Doppler 
velocity  LMB)  is  obtained  through 

VB)  ’  [£o  '  <±£b>]  I  <6-4) 

where  X  is  the  radio-wave  wavelength  2tt/$,  and  the  +  or  -  applies  to 
the  approaching  line  (first-order  band  near  f  )  or  receding  line,  re- 
spectively. 


B.  Statistical  Analysis  of  the  Doppler-Velocity  Estimator 

The  purpose  of  this  section  is  to  determine  whether  the  estimator 
is  unbiased  (its  expected  value  is  the  mean  current)  or  consistent  (its 
variance  about  the  mean  value  diminishes  as  the  coherent  integration  time 
T  increases) .  Only  the  approaching  line  is  analyzed  because  analysis 
of  the  receding  line  is  similar.  Combining  Eqs.  (6.3)  and  (6.4), 


or,  equivalently. 
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where  P^  is  the  power  spectrum  of  shifted  down  in  frequency  by  the 

amount  f  .  Extension  of  the  Parseval  theorem  (see  Appendix  C  for  the 
5 


proof)  yields 
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(6.5) 


where  the  *  denotes  the  complex  conjugate  and  the  angular  brackets  repre¬ 
sent  a  time  average  over  the  coherent  integration  time  T.  Here,  s(t) 

is  the  Bragg  signal  near  positive  f  multiplied  by  exp(-i2irf  t) .  From 

B  B 

Eq.  (5.2),  therefore, 

s(t)  =  B+  exp(i2frAt) 

where  A  is  the  current-induced  Doppler  frequency  and  B+  is  the  ampli¬ 
tude  of  the  Bragg-selected  wave  component  traveling  toward  the  radar.  As 
has  been  observed,  A  is  related  to  the  actual  depth-averaged  current 

A 

component  U (3)  through 


so  tnat 


s(t)  =  B  exp  [i4fr (U/X)  t] 


(6.6) 


The  subscript  of  B  has  been  dropped  for  simplicity. 

Equation  (6.6)  can  be  extended  to  the  condition  wherein  U  varies 
within  the  range  cell  probed  by  the  radar  by  separating  the  cell  into 
various  domains  within  each  of  which  0  is  approximately  constant.  The 
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linear  dimension  of  each  domain  must  be  large  compared  to  the  Bragg-se- 
lec'ed  wavelength  L  for  the  following  reasons. 


•  The  wave-number  selectivity  of  each  domain  must  remain 
sharp;  otherwise,  the  Doppler  spectrum  will  be  highly 
smeared — a  result  contradicted  by  observation. 

•  If  the  size  of  each  domain  is  smaller  or  comparable  to 
L,  the  horizontal  shear  of  the  current  cannot  be  ne¬ 
glected  in  the  hydrodynamics  and  dispersion  relationship 
derived  in  Chapter  II  is  not  valid. 


The  signal  s(t)  contains  contributions  from  all  the  domains, 


s(t) 
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a^B^  exp (i23u^t) 


(6.7) 


where  the  dependence  of  each  domain  on  area  a^  is  explicit.  Substitut¬ 
ing  this  model  of  s(t)  into  the  expression  for  the  centroid  estimator 
in  Sq.  (6.5)  yields 


(6.8) 


A 

Here,  A^  replaces  B^  exp(i2$U^t) .  Because  the  currents 
ent  domains  will  impart  a  large  differential  phase  to  the 
A.Aj*,  the  time-averaged  result  is  negligible  for  Z  /  k; 


in  the  differ- 
product  term 
therefore. 
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where 


A  straightforward  interpretation  of  the  current  estimate  is  now  the 
weighted  average  of  the  depth-averaged  current  field  over  the  range  cell 
and  over  the  coherent  integration  time  T.  To  make  this  interpretation 
more  explicit, 
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where 
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In  a  practical  context,  the  answer  to  the  question  of  an  unbiased  esti¬ 
mator  depends  on  the  desired  measurement  of  the  current  field;  it  is  ap¬ 
parent  that  this  measurement  should  correspond  to  the  current  averaged 
over  the  range  cell  and  over  T.  If  the  mean  wave-height  spectrum  is 
r.ot  uniform  over  the  cell,  the  above  weighted  average  will  differ  from 
the  desired  average  as  the  result  of  preferential  weighting  in  areas  of 

large  mean-squared  height.  If  the  spatial  variation  of  P  is  random 

m 

(uncorrelated  with  the  temporal  variation  of  P^)  then,  over  a  long  ob¬ 
servation  period,  will  tend  toward  the  desired  quantity. 

Reduction  in  the  random  variation  of  U  from  one  measurement  to 

c 

another  as  a  result  of  the  centroid  estimator  will  now  be  demonstrated. 
For  this  purpose,  spatial  averaging  can  be  ignored  which  will  greatly 
simplify  the  following  derivations;  therefore, 


1  fT/2 

U  =  -  dt  w  (t)  U (t) 

C  T  J- T/2 

where 
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w(t) 


P(t) 


P(t)  dt 


Assume  that  the  random  process  U(t)  is  statistically  independent 
of  the  normalized  spectral  process  w(t) .  Now,  w(t)  is  a  quotient  of 
two  random  variables;  the  denominator  is  a  time-averaged  version  of  the 
numerator.  Assuming  that  w  is  a  stationary  process,  then,  based  on  the 
following  definitions 

ElP(t)]  4  PQ 
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Var  [P  (t)  ]  “  o 


E{[p(t)  -  PQ][p(t  +  T)  -  PQ]|  k  Rp{T) 


and,  after  some  standard  manipulations. 


e[t  /p(t)  dt] '  po 


dtj  =  »p(T,(l-i) 


(6.12a) 


(6.12b) 


If  P  has  a  correlation  time  t  small  compared  to  T,  then 


Varj^  Jp(t)  dtj  =  20^  f 


(6.13) 


Because  the  power  spectrum  of  the  gaussian  process  (with  a  rectangular- 

window)  is  a  chi-square  of  two  degrees  of  freedom,  its  variance-to- 

squared-mean  ratio  is  unity.  The  variance  of  the  denominator  of  w(t), 

2 

therefore,  is  small  compared  to  P  and,  following  Barrick  [1978b], 


w(t)  = 


PQ(1  +  D) 
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(6.14) 


where 


After  some  lengthy  but  straightforward  manipulations.  an  expression 


for  the  variance  of  U  is 
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where  R^(t)  is  the  autocovariance  function  at  lag  t  of  the  random 

A  ■'*> 

process  U.  If  the  correlation  times  Ty,T^  f°r  t^ie  processes  U  and 
P,  respectively,  are  small  compared  to  T,  then 


Var [U  ] 
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when  t  >  t  .  If  t  «  t  ,  the  variance  of  U  is  approximately 
U  p  U  p  c 

doubled. 

The  above  calculations  demonstrate  the  variance  reduction  by  the 
centroid  estimator.  It  has  been  assumed  that  the  3ragg  width  is  caused 
by  the  statistical  character  of  the  current  field  and  wave-height  spec¬ 
tral  amplitude. 


C.  Experimental  Results  and  Interpretations 

It  has  been  demonstrated  that  the  centroid  estimator  in  Eq.  (6.3) 

/s 

can  be  used  to  produce  estimates  U(p)  of  the  depth-averaged  current  in 
the  direction  of  the  radar.  Eacn  U($)  is  actually  an  average  of  the 
current  component  over  the  coherent  integration  time  T  and  over  the 
entire  range  cell  probed  by  the  radar. 

A 

In  Figs.  6.1  through  6.4,  U  is  plotted  vs  time.  The  data  were 

collected  during  four  days  in  January  1978,  and  the  estimates  of  U  were 

obtained  using  the  approaching  Bragg  lines  because  the  receding  lines 

during  those  days  had  much  less  signal  strength.  Data  points  derived 

from  Bragg  lines  with  a  signal-to-noise  ratio  lower  than  6  dB  are  not 

plotted.  The  symbols  +  ,  x,  A,  and  >  denote  data  obtained  from  the 

four  range  cells  centered  at  9.75,  15.75,  21.75,  and  27.75  km,  respec- 

/\ 

tively.  Note  that  the  baselines  (U  =  0)  for  data  from  these  range 
cells  are  all  different;  each  has  been  displaced  upward  by  10  cm/sec 
for  each  increment  in  range,  as  indicated  by  the  solid  horizontal  lines. 

A 

On  January  19,  the  estimated  averaged  current  U  was  observed  to 
decrease  consistently  with  time  at  all  four  range  cells  and  radar  fre¬ 
quencies  during  uhe  four  hours  of  the  experiment.  Current  magnitude 
dropped  from  5  cra/sec  at  11; 27  local  time  to  zero  cm/sec  at  15:27;  su¬ 
perimposed  on  this  trend  was  a  small  statistical  variation  of  1  to  3  cm/sec. 
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d.  6.8  MHz 


Fig.  6.1.  RADAR- INFERRED  DEPTH -AVERAGED  CURRENT  AT  FOUR  RANGE  CELLS  VS 
TIME  DURING  19  JANUARY  1978.  Horizontal  lines  are  the  baselines  for 
the  different  range  cells. 
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Fig,  6.4.  RADAR-INFERRED  DEPTH-AVERAGED  CURRENT  AT  FOUR  RANGE  CfiLUS  VS 
TIME  DURING  25  JANUARY  1978 .  Horizontal  lines  are  the  baselines  for 
the  different  range  cells. 
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The  6.8  MHz  data  were  noisy  because  the  transmitted  pulse  was  reflected 
from  the  ionosphere  directly  overhead  and  experienced  Doppler  shifts 
caused  by  random  ionospheric  motions.  The  resulting  Doppler  spectrum 
was  relatively  broad  (a  few  tenths  of  a  hertz)  ana,  at  times,  the  Bragg 
lines  were  submerged  in  the  ionospheric  signals.  At  higher  radar  fre¬ 
quencies,  the  ionosphere  is  transparent  to  the  radio  waves. 

On  January  20,  the  radar-inferred  averaged-current  increased  stead¬ 
ily  with  time  at  all  range  cells  and  radar  frequencies;  the  total  increase 
was  approximately  5  cm/sec  over  the  four  hours  of  the  experiment.  Iono¬ 
spheric  interference  again  contaminated  the  data  at  6.8  MHz.  The  statis¬ 
tical  fluctuation  superimposed  on  the  temporal  trend  was  roughly  1  cm/sec. 

A 

On  January  24  and  25,  U  again  consistently  increased  with  time  at 
a] 1  range  cells  and  radar  frequencies  during  the  first  two  to  three  hours 
and  then  remained  constant  to  the  end  of  the  experiment;  the  total  in¬ 
crease  was  approximately  20  cm/sec  on  January  24  and  15  cm/sec  on  January 
25.  The  increase  of  current  with  radar  frequency  was  apparent  cn  those 
t^’O  days.  The  fluctuation  superimposed  on  the  temporal  trends  in  the 
estimated  current  was  again  roughly  1  cm/sec. 

The  observed  temporal  trend  in  the  depth-averaged  current  is  attrib¬ 
uted  to  two  sources — diurnal  and  semidiurnal  variations  in  the  tidal  cur¬ 
rent  component  [Sverdrup,  1942]  and  wind  speed.  On  January  19,  the  wind 
was  less  than  5  m/sec  and  came  from  the  west;  the  most  probable  explana¬ 
tion  for  the  steady  decrease  in  current,  therefore,  is  tidal  variation. 

On  January  24,  the  wind  rose  from  3  m/sec  at  11:20  local  time  to  6  m/sec 
at  15:40,  and  its  direction  shifted  from  030°  to  300°  T.  These  measure¬ 
ments  were  sporadic,  but  they  were  sufficient  to  indicate  a  rise  in  the 
wind  component  toward  the  radar.  The  increasing  wind  drag  on  the  ocean 
surface  may  have  been  partly  the  cause  of  the  increase  in  current  observed 
by  the  radar.  The  same  trend  occurred  on  January  25;  the  wind  rose  from 
4  m/sec  at  31:00  local  time  to  6  m/sec  at  14:28,  and  its  direction  shifted 
from  230°  to  300°  T. 

To  examine  the  correlation  between  tidal  currents  and  the  radar-mea¬ 
sured  current,  local  time  can  be  converted  into  a  phase  angle  from  the 
nearest  tidal  peak  as  recorded  in  the  tide  table  at  Point  Fort  in  San 
Francisco,  assuming  a  tidal  period  of  12  hours.  For  example,  the  data 
for  January  19  covered  phase  angles  from  approximately  100°  to  240°  and 
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the  data  for  January  24  covered  phase  angles  from  0°  to  125°.  It  is 
extremely  difficult  to  relate  water  movement  caused  by  tides  to  tidal 
current  [Lisitzin,  1974];  however,  the  nonoverlapping  phase-angle  ranges 
during  January  19  and  24  may  account  for  part  of  the  observed  reversal 
in  temporal  trends  in  the  currents. 
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Chapter  VII 


ESTIMATION  OF  DRIFT  PROFILE  FROM  CURRENT-INDUCED  DOPPLER  VELOCITIES 

As  has  been  discussed  in  the  previous  chapters,  the  estimated  cur¬ 
rent-induced  Doppler  velocity  is  a  weighted  average  of  the  ocean  current 
over  depth.  To  be  precise,  the  radar-inferred  Doppler  velocity  U(B)  at 
radio-wave  wave  number  B  has  been  shown  [Eq.  (6.1)]  to  be 

/■“ 

0(0)  =  46  J  Ug(z)  exp(-4Bz)  dz  (7.1) 

where  U^(z)  ls  the  component  of  the  horizontal  ocean  current  in  the 
direction  down  the  radar  beam  at  depth  z.  In  this  chapter,  the  depth 
variation  of  this  component  will  be  deduced  from  U(B)  measured  at  four 
radio-wave  wave  numbers  (6) • 

Theoretical  considerations  and  experimental  observations  [Bye,  1965; 
Shemdin,  1972;  Wu,  1975]  have  indicated  that  the  drift  profile  induced 
by  wind  drag  is  logarithmic, 


where 


U  =  surface  drift 
s 

U*  =  friction  velocity 
K  =  universal  Von  Karman's  constant 
zq  =  roughness  length 


(7.2) 


This  expression  is  not  valid  in  the  very  top  layer  of  the  ocean  where 
viscous  effects  cannot  be  neglected.  The  thickness  of  this  layer  is  on 
the  order  of  millimeters  [Bye,  1965].  The  importance  of  this  logarithmic 
profile  in  the  understanding  of  air/sea  interactions  leads  to  the  assump¬ 
tion  of  its  validity  in  a  first  attempt  to  invert  the  relationship  in  Eq. 
(7.1).  This  assumption  also  enables  determination  of  the  drift  profile 
from  in-situ  measurements  (see  Chapter  V)  so  that  comparisons  to  radar- 
inferred  current  can  be  made. 


83 


stag 


The  data  obtained  from  the  experiments,  however,  indicate  that  such 
assumptions  are  not  always  valid.  An  algorithm  is  derived,  therefore, 
for  inversion  of  the  Laplace  relation  in  Eq.  (7.1). 


A.  Comparisons  between  Radar-Inferred  Drift  Profile  and  In-Situ  Drift 
Measurements  for  Logarithmic  Profile 


The  logarithmic  profile  in  Eq.  (7.2)  is  assumed  to  be  valid  except 
for  a  very  thin  surface  layer  of  thickness  6  from  which  the  contribu¬ 
tion  to  the  integral  in  Eq.  (7.1)  is  negligible.  Introducing  a  variable 
s  =  4f3z, 


u(B) 


=  u 

s 


in 


exp (-436) 


£n  (s) 


exp(-s)  ds 


Because  6  is  much  smaller  than  the  radio  wavelength,  436  is  approxi¬ 
mately  zero  and  the  integral  (with  the  minus  sign)  in  the  above  equation 
is  then  a  constant;  actually,  it  is  equal  to  the  Euler  constant  y  =  0.5772 
[Carrier  et  al,  p.  186].  If  y  =  £n(r),  then 


0(3) 


(7.3) 


where  z  =  l/(4Br)  -  0.022  A,  and  A  is  the  radio-wave  wavelength.  In 
P 

other  words,  the  radar  probes  the  logarithmic  current  at  a  depth  of  ap¬ 
proximately  2  percent  of  the  radio-wave  wavelength.  If  a  linear  profile 
had  been  assumed  as  it  was  by  Stewart  and  Joy  [1974] ,  the  depth  probed 
by  the  radar  would  be  larger  than  z  by  a  factor  of  r  =  1.78. 

Based  on  this  result,  the  current-induced  Doppler-velocity  estimates 
0  at  four  radar  frequencies  discussed  in  Chaptsr  VI  can  be  interpreted 
as  the  actual  ocean  current  component  U  at  four  different  depths.  To 
compare  these  estimates  to  the  in-situ  measurements  obtained  by  drifters, 
the  radar-inferred  drifts  available  at  7-1/2  min  intervals  are  averaged 
over  a  time  period  covered  by  the  available  drifter  measurements.  The 
results  are  plotted  in  Figs.  7.1,  7.2,  7.3,  and  7.4  from  data  collected 
during  January  1978.  Note  that  these  plots  are  derived  from  positive 
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Fig.  7.1.  COMPARISON  OF  RADAR-INTERRED  DRIFT  AND  SPAR- 
BUOY  DRIFT  (CLUSTER  2)  OBTAINED  ON  19  JANUARY  1978. 
Horizontal  axis  is  logarithmic  in  depth,  and  a  loga¬ 
rithmic  profile  is  assumed.  The  approaching  Bragg 
line  was  used. 


Bragg  lines  (approaching  lines)  because  they  were  stronger  than  the  re¬ 
ceding  lines .  Superimposed  on  each  figure  are  the  currents  inferred  from 
in-situ  measurements  by  drifters  at  four  depths  of  immersion.  By  assuming 
a  logarithmic  current  profile,  the  measured  drifts  can  again  be  inter¬ 
preted  as  the  actual  current  at  four  different  depths  (see  Appendix  B) . 

From  these  plots,  the  following  consistencies  between  the  radar-in¬ 
ferred  drifts  and  the  in-situ  measurements  can  be  observed. 


The  same  change  of  current  with  time  (decreasing  the  first 
day  and  increasing  on  the  other  days)  is  detected  in  both 
measurements  at  the  various  depths  probed. 

The  decrease  of  current  with  depth,  in  most  cases,  is  seen 
in  both  measurements. 

On  January  25,  the  vertical  gradient  in  the  measured  current 
increases  with  time  from  both  the  radar  and  buoy  data.  This 
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Fig.  7.2.  COMPARISON  OF  RADAR- INFERRED  DRIFT  AND  SPAR- 
BUOY  DRIFT  (CLUSTER  1)  OBTAINED  ON  24  JANUARY  1978. 
Horizontal  axis  is  logarithmic  in  depth,  and  a  loga¬ 
rithmic  profile  is  assumed. 


can  be  explained  in  part  by  the  rise  in  wind  in  the  direc¬ 
tion  toward  the  radar,  as  described  in  Chapter  VI. 

It  is  also  apparent  from  these  plots  that  the  radar-inferred  current 
does  not  match  the  in-si tu  measurements  in  every  detail;  for  example,  the 
discrepancy  between  the  two  types  of  measurements  varies  from  less  than  1 
cm/sec  to  approximately  5  cm/sec.  Possible  reasons  for  this  discrepancy 
are  explained  below. 

1.  Resolution  Limitations 


The  resolution  of  the  radar  measurements  is  limited  by  the  fre¬ 
quency  resolution  of  the  Doppler  spectra.  As  discussed  in  Chapter  V,  the 


Fig.  7.3.  COMPARISON  OF  RADAR-INFERRED  DRIFT  AND  SPAR- 
BUOY  DRIFT  (CLUSTER  1)  OBTAINED  ON  25  JANUARY  1978. 
Horizontal  axis  is  logarithmic  in  depth,  and  a  loga¬ 
rithmic  profile  is  assumed. 


nominal  frequency  resolution  is  approximately  1  millihertz  which  results 
in  a  current  resolution  of  2.4,  1.2,  0.8,  and  0.6  cm/sec  for  radar  fre¬ 
quencies  at  6.8,  13.3,  21.7,  and  .°9.8  MHz,  respectively.  This  inability 
to  resolve  smaller  current  may  account  for  some  of  the  differences  ob¬ 
served  in  the  two  types  of  measurements. 

The  resolution  of  the  drifter  measurements  is  determined  by  the 
accuracy  of  the  range  measurements  which  are  a  few  meters.  Because,  over 
an  interval  of  an  hour,  this  range  uncertainty  results  in  a  velocity  un¬ 
certainty  very  much  less  than  1  cm/sec,  it  is  not  the  limiting  factor  in 
the  overall  accuracy  of  the  drifter  measurements. 

2.  Statistical  Uncertainties 

Because  of  the  temporal  and  spatial  averaging  inherent  in  the 
radar-inferred  current  estimator,  the  resulting  estimator  is  highly 
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Fig.  7.4.  COMPARISON  OF  RADAR- INFERRED  DRIFT  AND  SPAR- 
BUOY  DRIFT  (CLUSTER  1)  OBTAINED  ON  27  JANUARY  1978. 
Horizontal  axis  is  logarithmic  in  depth,  and  a  loga¬ 
rithmic  profile  is  assumed. 


consistent;  the  residual  statistical  fluctuation  after  the  temporal  trend 
is  removed  from  the  data  plotted  in  Figs.  6.1  through  6.3  is  on  the  order 
of  1  cm/sec.  The  drifter,  however,  senses  only  the  current  field  along  a 
specific  trajectory.  In  a  highly  variable  current  field,  drifter  measure¬ 
ments  can  differ  significantly  from  the  averaged  current  measured  by  the 
radar.  An  indication  of  current-field  variability  is  the  finite  width  of 
the  Bragg  lines  in  the  Doppler  spectra.  In-situ  evidence  of  this  varia¬ 
bility  is  apparent  in  the  discrepancy  between  the  motions  of  buoys  deployed 
6  km  apart  in  the  experiments.  The  current  profile  measured  by  the  two 
sets  of  spar  buoys  is  plotted  in  Fig.  7.5.  It  can  be  seen  that  the  cur¬ 
rent  magnitude  at  the  nearer  end  of  the  range  cell  is  smaller  than  that 
at  the  farther  end  by  at  least  5  cm/sec.  Larger  variability  in  drifter 
measurements  and  a  greater  discrepancy  between  radar-inferred  and  buoy- 
inferred  currents  have  been  reported  [Barrick  et  al,  1977b] . 
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Fig.  7.5.  VARIANCE  IN  DRIFTS  OF  SPAR  BUOYS  RELEASED 
IN  THE  SAME  RANGE  CELL.  Cluster  1  was  released  at 
25.8  km  from  Pescadero,  and  cluster  2  was  released 
at  20.5  km.  Data  were  collected  on  25  January  1978. 


This  current-field  variability  in  the  horizontal  plane  can 
generate  errors  in  estimating  the  vertical  gradient  (or  shear)  of  the 
current  from  the  in-situ  measurements.  Assume  that  a  wind-shear  flow 
produces  a  current  shear  that  is  uniform  over  the  range  cell.  Buoys  of 
varying  lengths  released  at  the  same  location  will  drift  simultaneously 
with  different  velocities  because  of  the  current  shear.  As  they  drift 
apart,  they  will  be  influenced  by  current  at  different  locations  and 
times.  Variations  in  the  drift  of  the  buoys  are  then  the  result  of 
current  variability  in  the  horizontal  plane  and  vertical  current  shear. 
This  may  explain  why  the  shear  observed  in  the  radar  data  differs  from 
that  in  the  drifter  data. 
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B.  Estimate  of  Friction  Velocity  U^ 

It  is  evident  from  Eq.  (7.2)  that  friction  velocity  U^  can  be 

estimated  from  the  slope  of  the  curve  U„  vs  £n(z).  The  radar  data 

P 

obtained  from  the  closest  range  cell  on  24  January  1978  are  used  here 
as  an  illustration  because,  as  can  be  seen  in  Fig.  7.6,  they  are  es¬ 
pecially  consistent  with  the  assumption  of  the  logarithmic  profile  [Eq. 
(7.2)].  The  slope  of  the  visual  fit  yields  =  1.6  cm/sec  based  on 
the  value  of  0.4  for  the  von  Karman  constant  K. 


Fig.  7.6.  LOGARITHMIC  CURRENT  SHEAR  MEASURED  BY  THE 
FOUR-FREQUENCY  RADAR.  Horizontal  axis  is  logarithmic 
in  depth,  and  the  logarithmic  profile  is  assumed  in 
the  inversion.  Vertical  bar  indicates  resolution 
limited  by  a  finite-time  Fourier  transform.  Data 
obtained  at  Pescadero  at  approximately  14:15  to  15: 

15  LT  on  24  January  1978;  range  6  to  13.5  km. 


During  this  particular  experiment,  the  wind  speed  V  down  the  radar 

beam  was  approximately  6  m/sec  which  resulted  in  a  wind  stress  of  (1/2) 

2  -3 

p  V  ,  where  p  is  the  air  density  (~1.2  x  io  gm/cc) .  The  drag  coef- 
a  a 

ficient  CD  defined  by 
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(7.4) 


Co  V 
Da 
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has  been  found  experimentally  to  be  approximately  1.3  x  10  ^  [Stewart, 
private  communication] ,  where  Pw  is  water  density  and  V  is  wind  speed 
at  an  anemometer  height  of  10  m.  The  anemometer  was  located  on  the  top 
of  the  mast  of  the  boat  used  to  deploy  and  track  the  buoys.  Its  height 
was  comparable  to  that  required  for  the  empirical  formula  [Eq.  (7.4)]. 

If  V  =  6  m/sec,  friction  velocity  is  -  0.8  cm/sec  which  is  one- 
half  that  measured  by  the  radar.  This  difference  may  be  caused  by  the 
uncertainties  in  the  wind  measurements  and  in  the  determination  of  the 
drag  coefficient. 
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C.  Drift-Profile  Estimation  from  Numerical  Inversion 

Inverting  the  integral  relation  in  Eq.  (7.1)  without  assuming  any 
specific  functional  form  of  current  vs  depth  distribution  is  complicated 
by  two  factors.  First,  the  value  of  the  current-induced  Doppler  velocity 
0(3)  is  known  only  at  four  values  of  the  argument  3  and,  second,  these 
known  values  are  inevitably  contaminated  with  measurement  noise.  The 
first  factor  constrains  the  number  of  quadrature  points  to  four  that  can 
be  used  to  approximate  the  integral  by  a  finite  sum.  To  minimize  the 
truncation  error  resulting  from  this  quadrature  approximation,  the  gaus- 
sian  quadrature  method  [Lanczos,  Chapter  VI;  Bellman  et  al,  Chapter  II] 
is  employed.  Both  truncation  error  and  measurement  noise  will  be  severely 
amplified  in  the  inversion,  and  this  always  occurs  when  inverting  the 
Fredholm  integral  equation  of  the  first  kind,  of  which  the  Laplace  rela¬ 
tion  here  is  an  example  [Phillips,  1962].  This  problem  is  solved  by  a 
stabilized  inversion  technique. 

1.  Reduction  to  Gaussian  Quadratures 

An  integral  can  be  approximated  by  summing  over  a  finite  number 
of  terms;  that  is. 


f  (xj 


w. 

X 


(7.5) 
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where  x^  are  the  quadrature  points  and  w^  are  the  weights.  In  the 
method  introduced  by  Gauss,  both  x^  and  w^  are  chosen  by  requiring 
that  Eq.  (7.5)  must  be  exact  for  any  polynomial  of  order  <2n-l. 

An  explicit  formula  for  w^  is  first  obtained  by  fitting  a 
polynomial  p^  ^ (x)  of  order  n  - 1  to  the  "data"  f (x)  at  all  n 
quadrature  points.  This  can  be  accomplished  by  forming  polynomials 
£L(x)  (i  =  1,2,..., n),  each  of  order  n-1  and  having  the  property 
that  it  is  zero  at  all  quadrature  points  except  x^;  that  is, 


Q. (x.)  =  6.  . 
i  D  ID 


(7.6) 


where  6..  is  the  Kronecker  delta.  Then, 
iD 


n 

?  .  (x)  =  X  f  (x. )  Q. 

n-1  +—L  a  i 

i=l 


(x) 


(7.7) 


will  match  f (x)  at  all  quadrature  points.  The  following  expression 
for  £L(x)  will  have  the  desirable  sampling  property  of  Eq.  (7.6): 


q  (x) 

Q.  (x)  =  — - 

i  q.(x.) 


(7.8a) 


where  q. (x)  =  (x  -  x, )  (x  -  x_) . . .  (x  -  x.  ,) (x  -  x. ,-)... (x  -  x  ).  Note 
l  12  l-l  l+l  n 

that  q.  (x)  is  a  polynomial  of  order  n-1  because  the  term  x-x. 

i  th  ^ 

is  missing  from  the  product  chain.  If  the  n  -order  polynomial  F^(x) 

is  defined  as 

F  (x)  =  (x-x.)  q.  (x) 
n  i  i 


an  equivalent  expression  for  (x) 


is 


.  F  (x) 
_  ,  ,  1  n 

2i(x)  F'(x.) 

n  i 


.  X-X. 

n  i  i 


(7.8b) 


which  is  the  Lagrange  interpolation  formula  that  yields  exact  values  of 
f (x)  at  any  point  x  if  f (x)  is  a  polynomial  of  order  n  -  1 .  Here , 
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the  prime  indicates  differentiation  with  respect  to  x,  and  F  (x) 

n 

n  zeros  located  at  the  n  quadrature  points  x^. 

Based  on  the  condition  that  Eq.  (7.5)  must  be  exact  for  p  Jx), 

n-i 

the  following  expression  for  the  gaussian  weights  w^  is  obtained: 


w.  =  I  Q.  (x)  dx  i  =  1,2,  _ _  n  (7.9) 

1  J-1  1 


This  equation  depends  on  the  exact  form  of  Q  (x)  which,  in  turn,  depends 
on  the  specification  of  the  quadrature  points  x^.  The  precise  values 
of  these  n  points  can  now  be  optimally  chosen  by  requiring  that  Eq. 
(7.5)  is  exact  for  all  polynomials  of  order  2n-l.  Equivalently,  Q  (x) 
should  be  such  that 


!  ' 


I 


w.  =  I  Q. (x)  dx  =  0  i  =  n  +  1,  . . . ,  2n 

1  -Li  1 

which  will  result  m  an  n-point  quadrature  formula  whose  accuracy  equals 
that  of  a  2n-point  quadrature  formula.  This  condition  is  analogous  to 
requiring  [Bellman,  1966,  Chapter  II]  that 


x1F  (x) 
n 


dx  =  0 


i  =  0,1, 


n  -  1 


and  this  orthogonality  condition  is  satisfied  by  the  Legendre  polynomial 

of  order  n.  The  n  quadrature  points  are  thus  identified  as  the  n 
fcl") 

zeros  of  the  n  -order  Legendre  polynomial.  With  x^  defined,  Q. (x) 
can  be  computed  and  used  to  determine  vr.  The  results  for  both  x^  and 
w^  are,  in  fact,  well  tabulated  [see,  for  example,  Abramowitz  et  al,  p. 
916]. 


To  reduce  the  integral  in  Eq.  (7.1)  to  quadratures,  the  trans¬ 
formation  x  =  2  exp(-4kQz)  -1  is  first  introduced  to  map  the  interval 
[0 ,°°]  into  [-1,1].  Here,  is  an  arbitrary  wave  number  chosen  to 

minimize  the  quadrature  error.  The  transformed  integral  can  then  be  re¬ 
duced  to  quadratures  as  described  above.  The  resulting  matrix  equation 
is 
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f  =  Ac 


where 


f  is  the  vector  wiuh  components 


(7.10a) 


2U (s.k  ) 

f.  =  - 3-2- 

1  s . 

i 


i  =  1,2, 3, 4 


(7.10b) 


c  is  the  vector  with  components 


=  Ug  |z (x J j  w_.  j  =  1,2, 3, 4  (7.10c) 


A  is  the  matrix  with  elements 


/II  \V 

Lij  (2  2  Xj) 


(7 . lOd) 


s.  is  the  normalized  wave  number 
1  1  0 

z(x.)  is  the  depth  defined  as 


z(x.)  =  - 


(i+  H) 


(7.10e) 


The  quadrature  points  x_.  and  weights  w_.  are  as  follows: 


0.8611363116 

0.3399810436 


0.3478548451 

0.6521451548 


Because  k^  is  chosen  as  the  smallest  of  the  four  wave  numbers  used  in 
the  experiments,  every  matrix  element  a^  is  bounded.  In  all  the  ex¬ 
periments  discussed  in  this  chapter,  k^  =  0.568  m  ^  which  corresponds 
to  the  radar  frequency  of  6.78  MHz.  From  Eq.  (7.10e),  therefore,  the 


depths  corresponding  to  the  gaussian  quadrature  points  x^  are  0.127, 
0.705,  1.952,  and  4.696  m. 

The  numerical  integration  indicated  in  Eq.  (7.10a)  was  pro¬ 
grammed  on  a  32-bit  computer  for  a  few  selected  drift  profiles,  and  the 
results  are  listed  in  Table  7.1.  It  can  be  seen  that  the  quadrature 
errors  are  no  more  than  3  percent  except  when  the  rapidly  varying  linear 
drift  profile  was  detected  to  a  relatively  great  depth  at  the  lowest 
radar  frequency. 

Table  7.1 

COMPARISON  OF  NUMERICAL- INTEGRATION  RESULTS  OBTAINED 
BY  GAUSSIAN  QUADRATURE  AND  EXACT  CALCULATIONS 


6  (m  1) 

U(0)  by 
Quadrature 

Exact 

A  - 

U($) 

Percent 

Error 

Uniform 

Current 

V2’  ■ 

20  cm/sec 

0.568 

20.00 

20.00 

0 

1.118 

20.00 

20.00 

0 

1.824 

20.00 

20.00 

0 

2.492 

20.00 

20.00 

0 

Linear 

Profile 

U 

e(z)  =  20  * 

(1  -  z)  cm/sec 

0.568 

-14.10 

-15.21 

7 

1.118 

2.05 

2.10 

3 

1 . 82-* 

9.04 

9.03 

0 

2.492 

11.98 

11.98 

0 

Exponential  Profile 
Up(z)  =  20  *  exp(-z)  cm/sec 


0.568 

1.118 

1.824 

2.492 


7.24 

10.56 

12.92 

14.27 


7.24 

10.56 

12.92 

14.27 


Logarithmic  Profile 


Uq (z)  =  20  +  2.5  *  £n(10z)  cm/sec 
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2. 


Direct  Inversion  and  the  Instability  Problem 

In  the  absence  of  truncation  errors  and  measurement  noise,  the 
drift  profile  Ug(z)  can  recovered  from  the  measured  Doppler  veloci- 
ties  U (p)  at  the  four  radio  wave  numbers  by  directly  inverting  Eq. 
(7.10a) ;  if  the  inverse  of  the  matrix  A  is  found  to  be  B,  then 


U 


6 


B.  .f  . 
13  3 


w . 

l 


i  =  1,  . . . ,  4 


(7.11) 


* 

i 
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For  the  matrix  inversion,  the  gaussian  elimination  method  (Dahlquist  and 
B^orck,  1974]  is  applied  and,  to  verify  the  algorithm,  the  numerical-in¬ 
tegration  results  in  Table  7.1  are  used  as  simulated  Doppler-velocity 
data  for  the  inversion.  The  results  of  the  inversion  are  exactly  the 
assumed  current  data  used  as  input  to  the  integration  routine,  as  ex¬ 
pected,  which  indicates  that  the  finite  register  length  of  the  computer 
does  net  introduce  significant  error. 

The  measured  Doppler  velocities  can  never  be  free  of  measure¬ 
ment  noise  £  which  may  be  caused  by  the  statistical  fluctuations  in¬ 
herent  in  the  random  nature  of  the  ocean  superimposed  on  noise  introduced 
by  the  radar  equipment.  The  quadrature  error  will  also  appear  as  a  noise 
component.  The  total  £  will  be  transmitted  to  the  inferred  profile 

UD(z)  through  the  matrix  inversion.  For  all  Fredholm  integral  equations 
p 

of  the  first  kind,  the  matrix  A  associated  with  the  kernel  is  unstable 
because,  in  f_  =  Ac_,  the  solution  vector  c_  is  very  sensitive  to  small 
perturbations  in  £  [Phillips,  1962].  All  measurement  noise,  therefore, 
will  be  highly  amplified  in  the  inversion  process,  and  the  inferred  pro¬ 
file  will  be  severely  distorted. 

To  demonstrate  the  instability  problem,  error  is  intentionally 
introduced  into  the  simulated  Doppler  velocities  in  the  following  manners 


f.  =  f. 
3  3 


|l  -  error  •  (-l)jJ 


where  the  error  is  0.01  in  the  first  case  and  0.1  in  the  second.  As  can 
be  seen  in  Table  7.2,  1  percent  noise  introduces  error  as  large  as  70 
percent  in  the  estimated  profile,  and  10  percent  noise  distorts  the  drift 
profile  beyond  recognition. 
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Table  7.2 


EXAMPLES  OF  INSTABILITY  IN  DIRECT  INVERSION.  Logarithmic 
profile  Ug(z)  =  20  -  2.5  *  £n(10z)  is  assumed. 


z 

(m) 

Actual  U„ 
(cm/sec) p 

Ug  from 

Direct  Inversion 

Noise 

Free 

1% 

Noise 

10% 

Noise 

0.13 

19.41 

19.41 

15.78 

-16.91 

0.71 

15.12 

15.12 

21.34 

77.36 

1.95 

12.57 

12.57 

3.66 

-76.50 

4. 70 

10.38 

10.38 

19.85 

105.16 

3.  Stabilized  Inversion 


The  matrix  equation  to  be  inverted  in  the  presence  of  measure¬ 
ment  noise  £  is 


f  =  Ac  +  e 


whose  least-square-error  solution  is  obtained  by  the  inversion 


t  t 

A  Ac  =  A  f 


where  t  indicates  that  the  matrix  A  is  transposed.  Because  A  is 
a  square  matrix  that  has  an  inverse,  this  least-square  solution  is  reduced 
to  solving  £  =  A£ — again  an  unstable  procedure. 

Further  restrictions  imposed  on  £  are  required  to  produce  a 
stable  solution.  Phillips  [1962]  was  the  first  to  propose  imposing  the 
condition  of  minimum  second  difference  of  £  in  solving  the  Fredholm 
integral  equation  of  the  first  kind.  Twomey  [1963]  simplified  this  ap¬ 
proach  and  suggested  minimizing  the  mean-squared  difference  between  c 


and  a  known  a  priori  profile  cq.  Bellman  et  al  [1966]  devised  algorithms 


for  the  above  constraints  but  from  a  dynamic -programming  standpoint  based 
on  an  iterative  approach  to  inverting  the  matrix  equations  considered  by 
Phillips  and  Twomey.  In  a  more  recent  survey,  Westwater  et  al  [1972]  de¬ 
scribed  a  statistical  technique  that  requires  as  inputs  the  covariance 
matrices  of  the  unknown  vector  c  and  measurement  noise  £. 
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In  the  method  employed  here,  a  specific  functional  form  (either 
logarithmic  or  linear)  is  first  assumed  for  the  current  vs  depth  distri- 


! 


bution,  and  a  first  estimate  of  the  current  at  four  depths  is  obtained. 
Based  on  a  least-square  fit,  the  current  is  then  extrapolated  at  the  four 
depths  dictated  by  the  numerical  quadrature  method.  This  extrapolated 
profile  is  then  converted  to  a  vector  cq  according  to  Eq.  (7.10c), 
which  is  then  used  to  stabilize  the  numerical  inversion  as  described 
below. 

The  formula  for  this  stabilized  inversion  can  be  derived  as 

follows.  A  solution  c_  is  sought  wherein  the  transformed  vector  Ac_ 

is  as  close  to  the  measurement  vector  f_  as  possible  and  the  deviation 

of  c  from  the  initial  solution  c.  is  small.  This  is  achieved  by 
—  -o 

minimizing  the  weighted  sum  of  the  mean-squared  errors, 


a '  ?  (?  V;  ‘  % 


+  X 


I 


(c.  -  c  .)' 
1  01 


Here,  the  weighting  factor  X  can  be  appropriately  called  a  Lagrange 
multiplier,  and  Q  is  a  quadrature  form  whose  minimum  is  attained  by 
differentiating  with  respect  to  any  component  c  of  c_  and  setting 
the  result  to  zero.  Therefore, 


?  (?  Vi '  £i)S  V**) +  x  l  <ei  *  v 


«  \k  -  0 


where  6^  is  the  Kronecker  delta  whose  sampling  property  can  be  used 
to  remove  the  summation  over  l.  As  a  result. 


S  2  4iaijcj  *  >ck  ”  S  <ifi  +  »=ok 


i  3 


In  matrix  form,  the  above  system  of  linear  equations  can  be 


written  as 


(AtA  +  Xl)  c  =  Afcf  +  Xcq 


(7.12) 
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Through  gaussian  elimination,  £  can  thus  be  obtained.  Note  that,  when 
X  vanishes,  Eq.  (7.12)  reduces  to  that  of  direct  inversion;  when  X  is 
exceedingly  large,  £  approaches  cq.  It  is  apparent  that  an  interme¬ 
diate  value  of  X  is  required  to  stabilize  the  inversion  to  improve  the 
initial  estimate  of  the  solution.  To  obtain  an  optimal  value  of  X, 
Phillips  [1962]  proposed  a  method  in  which  £  is  determined  from  Eq. 
(7.12)  and  the  root-mean-squared  error  e  =  (A£- f 5 t  (A£-  f_)  is  compared 
to  the  root -mean-squared  noise  e  in  the  measurement  vector  f .  The 
value  of  X  that  results  in  e  closest  to  the  noise  level  £  is  then 
the  optimal  value.  Phillips  reported,  however,  that  this  technique  was 
not  always  effective. 

In  the  approach  taken  here,  after  conversion  to  a  Doppler-ve- 
locity  unit  through  Eq.  (7.10b),  e  is  never  larger  than  1  cm/sec  in  both 
the  simulated  and  real  data;  it  increases  slowly  from  zero  at  X  =  0  to 
a  saturation  value  of  less  than  1  cm/sec  at  X  larger  than  unity.  The 
transition  at  X  =  0.001  is  fairly  sharp.  Although  the  inverted  vector 
solution  £  may  have  undergone  magnitude  changes  from  a  few  to  a  few 
hundred  cm/sec  as  X  drops  rrom  very  large  to  very  small  values,  the 
deduced  measurement  A£  varies  little  from  the  actual  measurement  f_. 

The  reason  for  this  behavior  is  the  unstable  nature  of  the  matrix  A.  As 
the  inverse  of  A  amplifies  measurement  noise,  the  matrix  itself  atten¬ 
uates  any  variations  in  £.  From  the  data  presented  in  Chapter  VI  and 
the  resolution  calculations  in  Section  A,  measurement  noise  here  is  also 
on  the  order  of  1  cm/sec.  Based  on  Phillip's  method,  a  value  of  X  =  0.1 
or  larger  is  required.  The  need  for  a  value  of  X  much  larger  than  unity 
implies  that  stabilized  inversion  will  not  significantly  improve  the  ini¬ 
tial  estimates  obtained  by  assuming  a  specific  functional  form  for  cur¬ 
rent  distribution  vs  depth. 

To  examine  in  detail  the  properties  of  the  inversion  algorithm, 
the  current-induced  Doppler  velocity  U($)  is  simulated  by  using  a  log¬ 
arithmic  or  linear  current  profile  Ug(z)  in  Eq.  (7.1)  for  which  the 
induced  Doppler  velocities  at  the  four  radio-wave  wave  numbers  are  the 
actual  current  at  known  depths,  as  discussed  in  Section  A.  A  fixed-per¬ 
centage  noise  can  be  added  to  the  Doppler-velocity  data  by  alternating 
the  sign  of  the  percentage  added  to  the  data  at  the  four  wave  numbers. 

To  generate  cq,  one  of  the  two  profiles  is  assumed  and  the  relationships 


in  Section  A  are  applied  to  relate  the  data  to  current  at  four  depths. 
Least-square  extrapolation  then  yields  cq  at  the  quadrature  points. 

In  Fig.  7.7,  the  inversion  algorithm  is  used  iteratively  to 
check  convergence  when  cq  is  intentionally  chosen  incorrectly  by  as¬ 
suming  the  wrong  profile.  In  Fig.  7.7a,  the  Doppler  data  contain  no 
noise  and  the  convergence  to  the  actual  current  profile  is  fairly  rapid. 
In  Fig.  7.7b,  5  percent  noise  was  added  to  the  data,  and  the  iterative 
approach  increased  the  deviation  of  the  derived  profile  from  the  actual 
profile.  It  is  apparent,  therefore,  that  the  best  approach  is  inversion 
without  iteration. 

Figure  7.8a  serves  as  a  check  for  the  algorithm.  Because  the 

data  are  free  of  noise  and  c  has  zero  weight,  the  inversion  is  essen- 

-o 

tially  perfect;  the  small  difference  between  the  resulting  current  pro¬ 
file  and  actual  profile  is  caused  by  propagation  of  the  quadrature  error 
which  must  have  been  even  smaller. 

In  Fig.  7.8b,  5  percent  noise  was  introduced  into  the  Doppler 
data.  Although  cq  has  the  correct  functional  form,  it  is  not  a  replica 
of  the  actual  current  because  of  the  added  noise.  The  inversion  output 
for  A  =  0.1  does  not  improve  the  extrapolated  profile,-  in  fact,  at  the 
last  quadrature  point  (4.7  m) ,  the  extrapolated  result  is  significantly 
closer  to  the  actual  profile  because  it  is  based  on  a  priori  knowledge 
of  the  functional  form  of  the  current  distribution. 

Figure  7.9  demonstrates  how  the  correct  choice  of  A  can  im¬ 
prove  the  extrapolation  results  based  on  the  assumption  of  a  wrong  pro¬ 
file.  This  is  encouraging  because  the  purpose  of  the  inversion  approach 
is  to  eliminate  the  need  for  a  priori  knowledge  of  the  functional  form 
of  current  distribution  vs  depth. 

In  Fig.  7.10,  the  inversion  method  is  applied  to  the  Doppler- 
velocity  data  derived  from  the  radar  backscatter  signal.  In  both  of  these 
figures,  current  is  again  assumed  initially  to  be  proportional  to  loga¬ 
rithmic  depth.  In  the  January  20  data,  the  logarithmic  assumption  is 
good,  and  the  inversion  result  is  highly  consistent  with  the  extrapolated 
profile.  In  the  January  24  data,  however,  there  is  a  sudden  depa\ture 
from  the  logarithmic  trend  at  the  longest  wavelength,  and  the  deduced 
current  profile  increases  slightly  with  depth  after  its  initial  sharp 
decrease. 
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b.  Five  percent  noise  added  and  Lag- 
range  multiplier  increased  to  0.1 


Fig.  7.8.  INVERSION  OF  SIMULATED  DOP¬ 
PLER  DATA,  LOGARITHMIC  PROFILE  ASSUMED 
IN  THE  INITIAL  ESTIMATES. 
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b.  Lagrangian  multiplier  A  =  0.01 


Fig.  7.9.  INVERSION  OF  SIMULATED  DOP¬ 
PLER  DATA  WITH  5  PERCENT  NOISE,  LOGA¬ 
RITHMIC  PROFILE  ASSUMED  IN  THE  INITIAL 
ESTIMATES . 
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Fig.  7.10.  INVERSION  OF  RADAR  DOPPLER' 
VELOCITY  DATA  AVERAGED  OVER  TIME. 


In  conclusion,  the  inversion  method  appears  to  work  reasonably 
well  given  the  fact  that  there  are  only  four  data  points.  If  a  signifi¬ 
cantly  larger  number  of  data  points  were  to  become  available  (a  radar 
capable  of  transmitting  and  receiving  sequentially  many  more  than  four 
frequencies) ,  the  inversion  will  be  more  reliable  because  the  random 
component  in  the  data  can  be  removed. 


Chapter  VIII 


BANDWIDTH  OF  THE  BRAGG  SIGNAL 


As  has  been  noted,  the  bandwidth  of  the  first-order  Bragg  signal  in 
the  radar  return  typically  ranges  in  units  of  current  from  5  to  15  cm/sec 
which  is  significantly  larger  than  the  resolution  limits  of  both  the  tem¬ 
poral  transform  and  spatial  filtering  that  result  from  the  backscattering 
process.  This  chapter  investigates  two  practical  methods  for  computing 
this  bandwidth  and  analyzes  its  behavior  as  such  radar  parameters  as 
pulse  width,  frequency,  and  time  of  sampling  and  other  uncontrolled  con¬ 
ditions  (such  as  wind  speed)  vary.  This  study  will  yield  valuable  insight 
into  the  ocean  processes  that  lead  to  Bragg-line  broadening. 


A.  Bandwidth  Estimators 

The  power  spectral  density  of  the  Bragg  signal  is  represented  by  P. 
at  f  =  i  •  Af,  where  Af  is  the  frequency  separation  between  adjacent 
spectral  points.  This  first-order  signal  is  assumed  to  be  nonzero  only 
in  the  frequency  region  determined  by  m^  <  i  <  m^;  all  summations  to 
be  discussed  in  this  chapter  will  cover  only  this  finite  range.  A  sec¬ 
ond-moment  bandwidth  w^  can  then  be  computed  as  follows: 


where 


(Hz) 


(8.1) 


S  = 

c  = 

A  = 


E  (i  -  c)  P.  =  second  moment  of  the  spectrum  about  the  centroid 


E  i  P./A  =  centroid  of  the  spectrum 
i  i 

E  P.  =  area  of  power  spectral  density  curve 
i  1 


This  estimated  width  is  the  frequency  separation  between  the  points  2.17 

dB  dov,n  from  the  peak  if  the  power  spectrum  is  exactly  gaussian  in  shape. 

The  major  weakness  in  this  estimate  is  its  sensitivity  to  noise  present 

at  the  frequency  bins  located  a  long  distance  away  from  the  signal  peak. 

2 

This  is  evident  in  the  weighting  factor  (i-c)  which  increases  as  the 
square  of  the  distance  from  the  centroid  c. 
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Area  bandwidth  is  a  much  less  sensitive  estimator, 

w  =  — —  •  Af  (Hz) 
max 


(8.2) 


where  A  is  the  area  of  the  spectral  curve  and  P  is  the  peak  value 

max 

of  the  power  spectrum.  This  estimate  is  the  exact  width  of  a  rectangular 
spectrum, 


(P 

I  max 


m,  <  1  <  nr 
1  -  ~  2 

otherwise 


The  principal  drawback  of  this  estimator  is  its  sensitivity  to  the  actual 
spectral  shape.  For  example,  when  a  band  of  signals  has  an  standing 
p  .ok,  w2  will  be  much  less  than  the  visual  width?  when  the  band  is 
fairly  uniform,  w2  will  be  approximately  equal  to  it. 

The  relation  between  w^  and  w2  can  be  derived  for  any  given 
spectrum.  For  example,  in  a  gaussian  curve. 


1.25 


and,  in  a  rectangular  spectrum. 


2.45 


In  the  gaussian  spectrum,  the  3  dB  width  is  larger  than  w^  by  a  factor 
of  1.18  and,  as  a  result,  it  is  well  approximated  by  w2- 


B.  Determination  of  First-Order  Region 

The  method  used  here  for  determining  the  cutoff  points  m  and  m 
that  define  the  first-order  Bragg  signal  differs  from  that  used  in  Chap¬ 
ter  V.  It  was  chosen  for  its  simplicity  and  because  it  yields  cutoff 
points  that  are  less  sensitive  to  spurious  signals  near  the  Bragg  lines. 

The  algorithm  on  which  this  method  is  based  is  as  follows. 


•  Obtain  a  smoothed  version  of  the  Doppler  spectrum  by  con¬ 
volution  with  a  rectangular  window  whose  width  is  10  fre¬ 
quency  bins. 

•  Search  for  the  peak  of  the  smoothed  spectrum  near  the  Bragg 
frequencies  ±ffi,  assuming  a  maximum  current  magnitude  of 
50  cm/sec . 

•  Locate  the  cutoff  points  10  dB  down  on  either  side  of  the 
peak. 


As  before,  these  cutoff  points  are  then  used  to  determine  the  desired 
parameters  from  the  unsmoothed  spectrum.  As  a  checx,  this  algorithm  was 
applied  to  the  data  set  discussed  in  Chapter  VI  to  obtain  the  radar-in¬ 
ferred  current,  and  the  resulting  current  estimates  differed  by  no  more 
than  1  cm/sec  from  the  estimated  currents  described  in  Chapter  VI.  Width 


estimates,  however,  are  more  sensitive  to  the  actual  locations  of  the 
cutoff  points. 

C.  Results 

The  set  of  radar  data  available  for  this  study  is  divided  into  the 
following  two  groups: 

•  data  collected  from  January  1975  through  January  1977 

•  data  collected  during  January  1978 


In  the  first  set,  the  four  radar  frequencies  employed  were  4.8,  6.8, 
13.3,  and  21.7  MHz  over  approximately  50  days  of  experiments.  During 
each  day,  four  1/2  hr  experiments  were  performed,  one  for  each  radar 
pulse  width  of  25,  50,  100,  and  200  ysec.  This  large  data  set  is  used 
to  study  the  behavior  of  the  Bragg  width  as  a  function  of  wind  speed  and 
direction,  radar  pulse  width,  distance  (range)  of  the  illuminated  ocean 
patch  from  the  radar,  and  the  radar  carrier  frequency,  respectively. 

The  second  set  has  the  same  data  base  used  in  Chapters  VI  and  VII. 
The  four  carrier  frequencies  were  6.8,  13.3,  21.7,  and  29.7  MHz  over 
five  days  of  experiments  (January  19,  20,  24,  25,  and  27).  Approximately 
four  hours  of  data  were  collected  each  day  at  a  fixed  radar  pulse  width 
of  50  ysec.  With  the  exception  of  January  20,  these  experiments  were 
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conducted  in  conjunction  with  simultaneous  in-situ  measurements  obtained 
by  spar  buoys  as  described  in  Chapter  V. 

This  second  smaller  data  set  is  used  to  determine  the  frequency  de¬ 
pendence  of  the  Bragg  width  over  a  slightly  different  range  of  frequen¬ 
cies,  and  the  temporal  trend  observed  in  the  ocean-current  estimate  (see 
Chapter  VI)  is  compared  to  the  temporal  behavior  of  the  Bragg  width  de¬ 
duced  from  this  data  set.  The  in-situ  measurements  of  the  total  current 
vector  are  also  employed  to  compute  the  Bragg  width  caused  by  the  finite 
radar  beamwidth,  and  compatibility  between  the  computed  and  observed 
widths  is  examined. 

The  estimated  width  is  denoted  by  w(a.,f.,L  ,r  ).  The  arguments 

i  j  k  m 

are  the  variables  defined  below. 


=  wind  speed  classified  in  four  categories 
i=l  from  0  to  4  knots 
i=2  from  5  to  9  knots 
i=3  from  10  to  14  knots 
i=4  higher  than  15  knots 


r 

m 


radar  carrier  frequencies,  increasing  with  index  j 

radar  pulse  widths,  increasing  with  index  k 

range  of  illuminated  ocean  patch  from  the  radar,  increasing 
with  index  m 


For  convenience,  the  estimated  width  w  in  units  of  frequency  is  con¬ 
verted  to  w  in  units  of  current  via 
c 

w  (cm/sec)  =  w  (Hz)  *  15000/f  (MHz) 
c  c 

The  wind  conditions  were  measured  at  the  Pescadero  radar  site  where 
an  anemometer  was  located  approximately  5  m  above  ground  level,  at  the 
S.E.  Farallon  Island,  and  at  Pigeon  Point  where  they  were  normally  re¬ 
ported  at  two-hour  intervals  by  the  weather  bureau  (Fig.  8.1).  The  wind 
speeds  at  these  three  locations  were  averaged  to  approximate  the  prevail¬ 
ing  winds  in  the  open  sea  covered  by  the  radar.  Variability  in  the  width 
estimates  was  reduced  by  averaging  over  those  obtained  on  different  days 
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Fig.  8.1.  LOCATIONS  OF  ANEMOMETERS.  These  are 
indicated  by  x  at  the  Farallons,  Pescadero,  and 
Pigeon  Point. 


but  under  the  same  classification;  the  resulting  width  is  denoted  by  w  . 
To  detect  trends  in  the  results,  further  analysis  was  based  on  the  as¬ 
sumption  that  the  variation  of  width  w  (a,f  ,L,r)  with  each  of  its  ar- 

c  c 

guments  was  separable  from  its  variation  with  the  other  arguments  in  the 
following  manner: 

wc<a,f,L,r)  =  «a(a)  wf(fc)  wL(L)  w^fr)  (8.3) 

Note  that  the  exact  functional  dependence  of  w  on  each  of  the  four  va¬ 
riables  is  not  specified  and,  as  will  be  demonstrated,  this  dependence 
on  most  of  the  variables  is  very  weak. 

By  normalizing  and  averaging,  assumption  (8.3)  enables  elimination 
of  the  dependence  of  wc  on  any  three  of  its  arguments  so  that  its  va¬ 
riation  with  the  remaining  argument  can  be  studied.  For  example,  the 
variation  of  Bragg  width  with  radar  frequency  f .  is 


W  (f.)  =  ~  ^  W  (f. 

c  j  64  .  4*  rkm  j 

•>  ■%  Lr  tn  “* 


(8.4) 


i,k,m 


where 


w.,  (f.)  = 

rkm  3 


w  (a.  ,f .  ,L  ,r  ) 
c  l  i  k  m 


w  (a. ,f  ,L  ,r  ) 
c  i  1  k  m 


Normalization  is  always  done  with  respect  to  the  width  at  the  index  value 
of  1  for  the  chosen  variable.  Associated  with  each  of  these  width  esti¬ 
mates  is  a  computed  standard  deviation  s;  for  example,  the  value  of  s 


at  f.  is 
3 


v  -  M. -  vvf} 


(8.5) 


In  Table  8.1,  all  the  width  estimates  in  the  first  column  are  equal 
to  unity  because  of  the  specific  normalizing  constant.  The  first  entry 
in  each  cell  is  the  normalized  area  width  derived  from  Eq.  (8.2);  the 
second  entry  is  the  second-moment  width  derived  from  Eq.  (8.1).  The 
bounds  of  each  estimate  are  established  by  the  standard  deviation  s. 
Note  that  both  the  area  and  second-moment  widths  generally  yield  similar 
values.  The  most  outstanding  exceptions  are  observed  in  the  row  for  the 
frequency-dependence  calculation  [third  row  for  wjfj];  these  discrep¬ 
ancies  are  presumably  caused  by  the  significant  fine-scale  differences 
in  the  Bragg-line  structures  at  the  different  carrier  frequencies. 

The  standard  deviations  are  large  compared  to  the  width  variations 
with  wind  speed,  range-cell  size,  and  range.  Conclusions  to  be  drawn  con¬ 
cerning  such  variations  are  determined  from  the  observed  average  behavior 
The  statistical  significance  of  these  conclusions  may  be  limited;  however 
width  variation  with  radar  frequency  exceeds  the  associated  standard  de¬ 
viation  which  indicates  a  conclusive  trend. 

The  behavior  of  the  Bragg  width  is  determined  from  both  data  sets, 
and  the  results  are  discussed  below. 
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Table  8.1 


NORMALIZED  WIDTHS  OF  BRAGG  LINES.  These  widths  are  derived  from 
the  radar  data  collected  from  May  1975  through  January  1977.  All 
data  with  less  than  a  20  dB  signal-to-noise  ratio  have  been  re¬ 
jected.  Eicner  the  receding  or  approaching  Bragg  line  was  used, 
depending  on  which  had  a  higher  signal  level.  The  first  entry  is 
normalized  area  width,  and  the  second  entry  is  normalized  second- 
moment  width.  The  bounds  are  set  by  the  standard  deviation  in 
each  estimate. 


i 

1 

2 

3 

4 

w  (a.)+ 

C  1 

1.00 

1.00 

1.03  ±  0.15 
1.03  ±  0.18 

0.99  ±  0.18 
0.98  ±  0.22 

1.05  ±  0.22 
1.11  ±  0.26 

w  (L. ) * 

C  1 

1.00 

1.00 

1.08  ±  0.18 
1.07  +  0.17 

1.04  ±  0.37 
1.05  ±  0.20 

1.03  ±  0.22 
1.04  ±  0.28 

w  (f.)* 

C  1 

1.00 

1.00 

C .91  ±  0.16 
0.97  ±  0.17 

0.59  ±  0.09 
0.76  ±  0.14 

0.45  ±  0.08 
0.62  ±  0.11 

A  ,  .  § 

w  (r . ) 

C  i 

1.00 

1.00 

1.01  ±  0.18 
1.02  ±  0.17 

1.07  ±  0.19 
1.07  ±  0.18 

1.10  ±  0.23 
1.13  ±  0.24 

^a.  =  wind-speed  categories  (0  to  4,  5  to  9,  10  to  14,  and  >15 
knots. 

*L.  =  radial  sizes  of  range  cell  (3.75,  7.5,  15,  and  30  km). 

*  1 

f.  =  radar  frequencies  (6.8,  13.3,  21.8,  and  29.8  MHz) . 

§  1 

r^  =  ranges  from  radar  (*=10,  15,  20,  and  25  km). 


3 .  Z  ipendence  on  wind  Conditions  (First  Data  Set) 

Although  the  last  wind-speed  category  (i  =  4)  was  set  aside 
for  all  velocities  higher  than  15  knots,  difficulty  in  maneuvering  the 
rhombic  antenna  at  speeds  higher  than  20  knots  set  an  upper  limit  on  the 
range  of  wind  conditions  observed.  The  results  indicate  that  the  width 
of  the  Bragg  line  increases  with  wind  speed  except  for  an  anomalous 
small  decrease  in  the  range  from  10  to  14  knots.  The  total  variation 
is  actually  only  approximately  10  percent  from  category  1  (i  =  1 )  to 
category  4  (i  =  4) . 

The  data  have  also  been  reclassified  according  to  wind  direc¬ 
tion,  irrespective  of  its  speed.  The  reference  direction  is  the  radar 


beam  axis  (approximately  315°  T) .  In 
Fig.  8.2,  a  cross  wind  is  indicated 
when  the  wind-velocity  vector  lies 
in  the  two  shaded  quadrants,  and  an 
up/down  wind  is  indicated  when  the 
vector  falls  in  the  unshaded  quadrants. 
Analysis  of  the  large  data  set  reveals 
that  the  Bragg  width  for  the  up/down 
wind  is  approximately  2  percent  larger 
than  that  for  the  cross  wind. 

2.  Dependence  on  Radar  Pulse  Width  (First  Data  Set) 

For  the  radar  pulse  widths  used  in  the  experiments,  the  corre¬ 
sponding  radial  sizes  of  the  illuminated  ocean  patch  are  3.75,  7.5,  15, 
and  30  km.  The  results  in  Table  8.1  show  that  the  width  of  the  first- 
order  line  increases  by  8  percent  when  the  radial  extent  of  the  illumi¬ 
nated  area  increases  from  3.75  to  7.5  km;  a  further  increase  apparently 
causes  a  small  reduction  in  the  Bragg  width. 

3.  Dependence  on  Range  from  Radar  (First  Data  Set) 

The  four  range  cells  sampled  by  the  receiver  are  approximately 
10,  15,  20,  and  25  km  away  from  the  radar.  Analyses  indicate  that  the 
Bragg  width  increases  gradually  with  range. 

4.  Dependence  on  Radar  Frequency 
a.  First  Data  Set 

Both  area  and  second-moment  widths  decrease  sharply  with 
increasing  radar  frequency  although  the  second-moment  estimate  indicates 
a  relatively  weaker  frequency  dependence.  Note  that  the  width  estimates 
prior  to  normalization  are  in  units  of  cm/sec.  If  the  frequency  depen¬ 
dence  of  the  Bragg  width  is  approximated  bv  a  power  law  (w  ~  f  a) ,  the 

c  c 

data  in  Table  8.1  can  be  used  to  estimate  a.  Based  on  the  tabulated 


Fig.  8.2.  WIND-DIRECTION 


results,  a  =  0.5  for  area  width  and  0.3  for  the  second-moment  width. 


b. 


Second  Data  Set 


j 

i 


The  nominal  values  of  the  radar  carrier  frequencies  used 
in  these  experiments  were  6.8,  13.3,  21.7,  and  29.75  MHz.  After  analyz¬ 
ing  the  data  from  the  first  range  cell,  an  approximate  power-law  depen- 
-ct 

dence  (w  ~  f  )  was  observed  with  a  -  0.4. 
c  c 

5.  Correlation  with  Radar- Inferred  Current  (Second  Data  Set) 

Although  the  radar-inferred  current  component  examined  in  Chap¬ 
ter  VI  showed  a  consistent  trend  with  time  during  each  of  the  five  days 
of  experiments,  the  Bragg  width  did  not  exhibit  the  same  behavior.  For 
example,  on  January  25,  the  estimated  radial  current  increased  from  0  to 
20  cm/sec  during  the  four  hours  of  experiments,  and  the  Bragg  width  [es¬ 
timated  using  Eq.  (8.2)]  fluctuated  about  the  mean  value  of  5.7  cm/sec 
at  the  radar  frequency  of  13.3  MHz. 


6.  Correlation  with  Width  Estimated  from  In-Situ  Current 


in  the  area  illuminated  by  the  radar  beam  was  not  uniform  (see  Fig.  7.5); 
the  turbulent  nature  of  these  currents  has  also  been  reported  by  Sverdrup 
[1942] .  Assuming  that  the  ocean  current  within  a  given  range  cell  is 
turbulent  and  that  a  standard  deviation  0 is  associated  with  its  prob¬ 
ability  distribution  function,  the  width  of  the  Bragg  line  is  also  ap¬ 
proximately  o^.  The  in-situ  measurements  did  not  provide  sufficient 
data  from  which  a  statistically  significant  estimate  of  0  could  be 
extracted;  however,  the  difference  in  the  measurements  by  the  two  sets 
of  buoys  was  of  the  same  order  of  magnitude  as  the  width  of  the  Bragg 
line. 


{ 


D.  Summary 

This  section  summarizes  the  characteristics  of  the  Bragg-line  width 
(in  units  of  cm/sec)  as  obtained  from  processing  the  HF  ocean-backscatter 
radar  data  consisting  of  nearly  10,000  Doppler  spectra. 


•  It  is  inversely  proportional  to  radar  carrier  frequency  f  ; 
a  power-law  fit  to  the  data  yielded  a  frequency  dependence 
of  «f_0-4. 
c 


i 

i 


t 


< 


•  It  is  a  weak  function  of  the  prevailing  wind  speed  in  the  J 

range  from  0  to  20  knots.  Only  up  to  a  10  percent  increase  i 

was  observed  when  the  wind  speed  rose  at  least  threefold.  j 

•  when  wind  blows  across  the  radar  beam,  the  Bragg  width  is  [ 

approximately  2  percent  larger  than  when  wind  blows  along  j 

the  beam.  I 

•  It  increases  by  8  percent  when  the  radial  extent  of  the  il-  | 

luminated  ocean  patch  increases  from  3.75  to  7.5  km;  further  | 

increases  to  15  and  then  to  30  km  apparently  will  reduce  the  I 

Bragg  width  by  a  few  percent.  J 


It  increases  by  approximately  10  percent  as  the  center  of 
the  range  cell  sampled  by  the  radar  recedes  from  10  to  25 
km  as  measured  from  the  radar.  This  broadening  is  most 
likely  caused  by  a  reduction  in  the  signal-to-noise  ratio 
with  range  which  results  in  a  higher  probability  of  noise 
corruption  in  the  width  estimator. 


It  is  uncorrelated  with  the  radial  component  of  the  ocean 
current  detected  by  the  radar. 


117 


It 

t-s 


Any  attempt  to  explain  the  above  observed  features  of  the  Bragg 
signal  will  require  a  model  of  the  underlying  processes  that  cause  the 
Bragg  width.  The  possibility  that  broadening  is  the  result  of  radar 
system  artifacts  (such  as  antenna  motion)  has  been  considered  in  Chapter 
V;  however,  experimental  evidence  indicated  that  very  little  broadening 
could  be  accounted  for  in  this  manner.  The  effect  of  a  finite  beamwidth 
illuminating  a  uniform  current  has  also  been  examined,  and  experimental 
evidence  revealed  that  this  accounted  for  part  of  the  observed  width  and 
also  part  of  the  width  reduction  with  frequency. 

In-situ  measurements  showed  that  the  waters  off  the  San  Francisco 
peninsula  have  an  eddy-current  component  which  is  random  in  nature.  If 
this  component  does  broaden  the  Bragg  line,  the  following  characteristics 
must  be  inherent  in  the  eddy-current  structure. 


•  The  radial  correlation  length  of  the  random  current  field 
must  be  much  less  than  3  km  to  account  for  the  small  change 
in  the  Bragg  width  with  the  radial  extent  of  the  observed 
ocean  patch. 

•  The  dependence  of  the  Bragg  width  on  carrier  frequency  in¬ 
dicates  one  or  both  of  the  following  possibilities. 

-  The  turbulent  current  field  has  a  standard  deviation 
that  increases  with  depth.  As  radar  frequency  drops, 
the  effective  depth  probed  by  the  radar  increases  and 
the  Bragg  line  broadens. 

-  The  azimuthal  correlation  length  of  the  random  current 
field  is  comparable  to  the  azimuthal  size  of  the  illum¬ 
inated  ocean  surface.  At  15  km,  the  azimuthal  size 
ranges  from  3  to  13  km  for  the  radar  frequencies  used. 

•  The  weak  correlation  of  the  Bragg  width  with  wind  speed 
suggests  that  the  principal  mechanism  for  the  generation 

of  current  turbulence  is  not  interactions  with  the  prevail¬ 
ing  wind. 


Further  experimental  work  is  required  to  validate  the  above  observations. 
Theoretical  studies  of  the  effect  of  small-scale  horizontal  variations 
of  ocean  current  on  wave  phase  velocity  may  also  clarify  the  broadening 
phenomenon  of  the  Bragg  line. 


Chapter  IX 


CONCLUSIONS  AND  RECOMMENDATIONS 


A.  Conclusions 

The  feasibility  of  remotely  sensing  ocean  surface  current  and  its 
distribution  with  depth  using  a  four-frequency  HF  backscatter  radar  has 
been  studied  theoretically  and  experimentally.  Theoretical  investiga¬ 
tions  indicated  that  the  result  obtained  by  Stewart  and  Joy  [1974]  for 
ocean-wave  phase  velocity  in  the  presence  of  surface  current  with  a  non¬ 
zero  vertical  gradient  is  valid  if  the  vertical  gradient  is  small.  Their 
result  was  extended  to  second  order  for  an  exponential  current  profile; 
however,  for  conditions  prevailing  during  the  experiments  reported  here, 
this  second-order  effect  was  not  important.  A  Bragg-scattering  theory 
based  on  simple  wave -propagation  concepts  was  derived,  and  it  was  demon¬ 
strated  that  the  current  component  in  the  direction  of  the  radar  can  be 
estimated  from  the  position  of  the  Bragg  signal  in  the  ocean  Doppler 
spectrum. 

Data  collected  by  the  Stanford  four-frequency  radar  located  at  Pes- 
cadero  on  the  California  coast  were  analyzed  to  generate  an  estimate  of 
the  ocean  current  component  in  the  direction  of  the  radar.  A  15  min  co¬ 
herent  integration  time  produced  resolutions  of  2.4,  1.2,  0.8,  and  0.6 
cm/sec  at  the  radar  frequencies  of  6.8,  13.3,  21.7,  and  29.8  MHz.  Exper¬ 
iments  during  January  1978  revealed  that  fluctuations  in  the  current 
estimates  were  generally  no  more  than  1  to  3  cm/sec  superimposed  on  steady 
temporal  trends. 

Spar  buoys  1,  3,  6,  and  12  m  long  were  deployed  and  tracked  by  a 
microwave  system  to  obtain  in-situ  measurements  of  the  ocean  surface  cur¬ 
rent  which  were  then  compared  to  the  radar-inferred  current  component  by 
assuming  that  it  varies  in  proportion  to  logarithmic  depth.  The  discrep¬ 
ancy  between  the  two  types  of  measurements  was  no  larger  than  the  dis¬ 
crepancy  between  the  in-situ  measurements  made  along  two  different  tra¬ 
jectories  within  the  same  range  cell  probed  by  the  radar.  The  general 
behavior  of  the  ocean  current  with  time  and  depth  measured  by  the  spar 
buoys  agreed  with  that  measured  by  the  radar. 


Reconstruction  of  the  ocean-current  distribution  with  depth  from  the 
radar  Doppler  measurements  without  assuming  an  explicit  functional  form 
of  the  current  profile  requires  inversion  of  a  Fredholm  integral  equation 
of  the  first  kind.  A  numerical  inversion  technique  was  investigated, 
using  a  stabilization  scheme  based  on  an  initial  estimate  of  the  current 
profile.  The  implemented  algorithm  indicated  improvement  of  this  esti¬ 
mate  when  applied  to  simulated  Doppler  data. 

From  a  large  set  of  radar  data  collected  from  May  1975  through  Jan¬ 
uary  1978,  several  characteristics  of  the  Bragg  width  were  analyzed.  This 

width  (in  units  of  current)  was  observed  to  vary  with  radar  frequency  f 

-a 

as  f  (with  a  -  0.3  to  0.5),  and  its  dependence  on  radar  pulsewidth, 
wind  speed,  or  wind  direction  was  much  weaker — up  to  a  few  percent  change 
for  the  various  conditions  encountered. 

This  research  has  shown  that  HF  backscatter  with  average  RF  power  of 
a  few  watts  can  probe  the  radial  component  of  horizontal  current  and  its 
distribution  with  depth  in  an  ocean  patch  tens  of  kilometers  away.  Reso¬ 
lutions  of  a  few  cm/sec  can  be  obtained  during  time  intervals  of  a  few 
minutes  in  an  ocean  patch  of  linear  dimension  on  the  order  of  kilometers. 
The  current  component  orthogonal  to  the  radar  beam  can  be  monitored  by  a 
similar  radar  deployed  at  a  different  location.  The  ocean  surface  can  be 
scanned  by  range  gating  and  beam  steering  to  achieve  simultaneous  coverage 
of  large  areas.  These  features  of  the  radio  technique  are  very  attractive 
when  compared  to  conventional  oceanographic  measurements. 

B.  Recommendat ions 

The  radar  system  considered  here  employed  a  fixed  beam.  To  map  ocean 
current  within  a  full  180°  field  of  view,  the  radar  has  been  modified  to 
receive  the  signal  on  an  array  of  eight  loop  antennas,  and  data  have  been 
collected  by  this  new  system  [Teague  et  al,  1978].  A  challenging  area  of 
research  is  in  devising  an  optimal  scheme  for  combining  the  signals  on 
the  eight  antennas  so  that  the  antenna  beam  can  be  steered.  Beam  forming 
by  conventional  spatial  Fourier  transform  and  by  maximum  entropy  [Burg, 
1975]  should  be  compared.  The  method  first  suggested  by  Crombie  [1972] 
and  extended  by  Barrick  [1976]  to  determine  the  angle  of  arrival  by  com¬ 
paring  the  phases  of  signals  from  the  various  receiving-array  antennas  is 
also  applicable  here. 


To  further  understand  air/sea  interactions,  it  is  desirable  to 


measure  the  distribution  of  ocean  current  with  depth,  especially  when 
the  prevailing  wind  speed  is  high  (over  20  knots) ,  Such  data  did  not 
become  available  partly  because  of  the  prevailing  weather  conditions  on 
the  California  coast  during  the  experiments  and  partly  because  of  the 
difficulty  in  handling  the  rhombic  transmit  antenna  during  high  winds. 
The  new  system  described  by  Teague  [1978]  replaces  the  balloon-supported 
rhombic  antenna  with  vertical  whips  supported  by  a  fixed  mast.  The  radar 
system  can  now  operate  in  high  winds  and  obtain  current-shear  measure¬ 
ments  as  a  function  of  wind  conditions. 

In  this  work,  it  was  assumed  that  the  variation  of  ocean  current  in 
the  horizontal  plane  could  be  ignored.  This  assumption  is  violated  when 
the  scale  of  variation  is  comparable  to  or  smaller  than  the  wavelength 
of  the  ocean  waves.  To  incorporate  the  effects  of  such  horizontal  cur¬ 
rent  structures  into  the  dispersion  relation,  tedious  derivations  are 
required;  however,  the  results  may  account  for  part  of  the  Bragg-line 
broadening  phenomenon. 


Appendix  A 


DEGREE  OF  FREEDOM  OF  HAMMING-WINDOWED  SPECTRAL  ESTIMATOR 


Consider  a  real-time  series  s(t)  of  duration  T.  Assume  that 
s (t)  is  a  zero  mean  stationary  process  with  an  autocovariance  function 
Y(u)  defined  as 

Y(u)  =  E  [s (t)  s (t  +  u)] 

Denote  the  time  window  by  q(t).  This  function  is  zero  outside  the  in¬ 
terval  [-T/2,T/2].  The  windowed  power  spectral  density  estimator  P(f) 
at  frequency  f  is 

1  I  r  2 

P(f)  =  —  I  s(t)  a  (t)  exp(-i2irft)  dt 


The  integration  limits  [-°°,“]  are  assumed.  Taking  the  ensemble  aver¬ 
age, 


E  [P  (f )  ]  =  ^  JJ  E  [s  (t)  s(t')]  q  (t)  q(t')  exp  [-i2irf  (t  -  t')]  dt  dt' 
=  j  Y(u)  w(u)  exp(-i2irfu)  du 


(A.  1) 


where 


w(u)  =  ^  J 


q(t)  q (t  -  u)  dt 


The  above  manipulations  establish  the  relationship  between  the  lag 
window  w(u)  used  by  Jenkins  and  Watts  [1968]  and  the  time  window  q(t). 
By  comparing  Eq.  (A.l)  to  their  equations  (6.3.28)  and  (6.3.29),  it  can 
be  seen  that  w(u)  is  simply  w(u)  normalized  by  w(0).  The  equivalent 
degree  of  freedom  V  can  be  computed,  therefore,  by  using  their  equation 
(6.4.17) , 


where 


The  Hamming  window  is  a  complete  cycle  of  a  cosine  waveform  on  an  8 
percent  pedestal, 

(0.54  +  0.46  cos  (2irt/T)  |t|  <  T/2 
0  otherwise 

The  equivalent  lag  window  w(u)  can  thus  be  calculated  by  numerical  in¬ 
tegration,  and  the  degree  of  freedom  for  the  Hamming  window  is  approxi¬ 
mately  4. 


Appendix  B 


SPAR-BUOY  DRIFT  IN  A  CURRENT  THAT  VARIES  WITH  DEPTH 

The  following  analysis  is  based  on  the  work  of  Bye  [1965] .  It  is 
included  here  for  completeness. 

Let  SL  be  the  depth  of  immersion  of  the  spar  buoy  and  U(z)  be  the 
horizontal  current  that  varies  with  depth  z  below  the  ocean  surface. 
Assume  that  the  buoy  will  drift  along  with  the  current  at  a  speed  denoted 
by  UD. 

For  the  type  of  flow  here,  the  Reynold  number  is  on  the  order  of  a 

few  thousand,  and  the  drag  coefficient  is  independent  of  the  velocities 

involved  over  this  range  [Kaufmann,  1963].  The  drag  force  acting  on  any 

2 

part  of  the  spar  buoy  is  thus  proportional  to  [UD~U(z)]  ,  with  a  sign 
that  follows  the  sign  of  [U^-UC^.)]. 

For  the  buoy  to  remain  in  a  vertical  position  (stable  equilibrium) , 
the  drag  force  should  distribute  itself  so  that  its  effect  on  the  buoy 
is  zero.  In  this  event,  there  must  be  a  critical  depth  z^  within  the 
interval  [0,£]  at  which  the  drag  force  is  zero;  in  other  words, 

U(z  )  =  U^  (B.l) 

C  D 

The  requirement  for  equilibrium  can  thus  be  written  as 


(B.2) 


whose  solutions  are  sought  for  both  the  linear  and  logarithmic  profiles. 


If  the  profile  is  linear  [U(z)  =  UQ  +  az],  direct  substitution  into 
Eq.  (B.2)  yields 


The  obvious  solution  is  z^  =  J2./2.  (The  other  two  solutions  are  imagi¬ 
nary.)  For  a  linear  current  profile,  therefore,  the  buoy  drifts  with  a 
velocity  equal  to  the  current  at  a  depth  one-half  of  its  immersed  length. 

If  the  profile  is  logarithmic  [U(z)  =  Uq +  •  £n(z)],  direct  sub¬ 

stitution  into  Eq.  (B.2)  results  in 


which,  if  r 


dz 


z/z  and  r„ 
c  0 


1/7.  , 
C 


becomes 


£n2  (r)  dr 


£n2(r)  dr 


Two  integrations  by  parts  obtain 


“=  *n2(r0)  [*n(V  “  2]  +  2 

Solving  by  successive  iterations,  =  0.2725  or  z^  =  0.2725  1.  For 

a  logarithmic  current  profi?.e,  therefore,  the  buoy  drifts  with  a  velocity 
equal  to  the  current  at  a  depth  slightly  over  one-fourth  of  its  immersed 
length. 
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Appendix  C 


PARSEVAL'S  THEOREM  AND  ITS  EXTENSIONS 


Given  a  signal  s<t)  in  the  time  interval  [-T/2,T/2],  a  conven¬ 
tional  assumption  used  in  deriving  its  discrete  Fourier  transform  (DFT) 
is  that  it  repeats  itself  outside  this  interval.  Based  on  this  assump¬ 
tion,  the  signal  can  be  written  as  a  Fourier  series  with  a  fundamental 
frequency  Af  =  T  7  The  sampling  frequency  f^  is  assumed  to  be  larger 
than  twice  the  bandwidth  of  s(t).  Within  the  Nyquist  frequency  interval 

of  [-f  /2,f  /2],  therefore,  the  Fourier  series  coefficients  a  are 
s  s  n 

the  same  as  che  DFT  coefficients  S(nA^),  and  the  summations  below  are 
restricted  to  this  interval.  As  a  result, 

s  (t)  =  ^  S(nAf)  e:cp(i2imAft) 
n 

and 

S(nAf)  =  ^  f s(t)  exp(-i27rnAft)  dt 


where  the  integration  limits  are  assumed  to  be  from  -T/2  to  T/2. 

Parseval's  theorem  relates  the  energy  in  s(t)  to  that  in  the  DFT 
coefficients.  The  average  energy  in  s  (t)  is 


■7*11 

^  m,n 


S(mAf) 


S*(nAf) 


exp|i27r(m  -  n) 


where  *  denotes  the  complex  conjugate.  Because  the  complex  exponentials 
exp (i2mn £^t;  (m  is  an  integer)  are  orthogonal  in  [-T/2,T/2],  Parse- 
val ' s  theorem 

T  1  j  | s (t )  j2  dt  =  ^  P(nAf) 

^  n 

2 

follows,  where  P(nA^)  =  |s(nA^)|  is  the  power  spectral  density. 

A  second  theorem  relates  the  first  moment  of  the  power  spectrum  tc 
an  integral  in  the  time  domain.  If  differentiation  with  respect  to  time 
t  is  denoted  by  a  prime,  then 


IT  I  s (t)  s* ' (t)  dt 


dt  2_,2-,  27rn^fS^f)  S*(nAf) 
m,n 


exp^i27T(m  -  n)  AftJ 


by  direct  manipulation  of  the  Fourier  expansion  of  s(t).  Again  using 
the  orthogonal  relationship  for  the  complex  exponentials, 

mA^P(mA^)  =  (2ttT)  ^  J"  is(t)  s*'(t)  dt 


A  third  theorem  relates  the  second  moment  of  the  power  spectrum  to 
an  integral  in  the  time  domain.  The  following  relcition  can  be  obtained 
from  the  above  type  of  manipulations: 

^  (nAf)2  P(nAf)  =  (2t r)  2T_1  J  s' (t)  s*'(t)  dt 


m$0mm 


Appendix  D 


RELATIONSHIP  BETWEEN  THE  WINDOWED  TRANSFORMS 
OF  A  FUNCTION  AND  ITS  DERIVATIVE 


Let  the  function  be 
of  its  derivative  f ' (x) 


f  (x)  and  define  the  windowed  Fourier  transform 
by 


F1;s) 


f ' (x)  exp(-isx)  dx 


Here,  the  window  is  assumed  to  be  zero  outside  the  interval  [-L/2,L/2]. 
The  integration  interval  is  assumed.  Integrating  by  parts. 


I  +00 

F- (s)  =  w{x)  f (x)  exp(-isx) 

J_  —  00 


+  is  I  w{x)  f(x)  exp(-isx)  dx 


is  j  w 
"/  W 


'  (x-  f (x!  exp(-isx)  dx 


The  first  term  is  zero  because  the  window  function  w(x)  is  nonzero  only 
in  [-L/2,L/2].  The  third  term  is  negligible  compared  to  the  second  be¬ 
cause  the  width  of  w(x)  is  large  compared  to  the  wavelength  2tt/s  of 
the  complex  sinusoid.  As  a  result. 


F1(s)  =  isF(s) 


where  F(s)  is  the  windowed  Fourier  transform  of  f (x) . 
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